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Abstract 

The numerical studies of (ij.2_j^2-wave pairing in the two-dimensional (2D) and the 2-leg 
Hubbard models are reviewed. For this purpose, the results obtained from the determinantal 
Quantum Monte Carlo and the density-matrix renormalization-group calculations are presented. 
These are calculations which were motivated by the discovery of the high-Tj, cuprates. In this 
review, the emphasis is placed on the microscopic many-body processes which are responsible for 
the c?2,2_j^2-wave pairing correlations observed in the 2D and the 2-leg Hubbard models. In order 
to gain insight into these processes, the results on the effective pairing interaction as well as the 
magnetic, density and the single-particle excitations will be reviewed. In addition, comparisons 
will be made with the other numerical approaches to the Hubbard model and the numerical 
results on the t-J model. The results reviewed here indicate that an effective pairing interaction 
which is repulsive at (tt, tt) momentum transfer and enhanced single-particle spectral weight 
near the (7r,0) and (0,7r) points of the Brillouin zone create optimum conditions for d^2_y2- 
wave pairing. These are two effects which act to enhance the d^^-yi-wave pairing correlations 
in the Hubbard model. Finding additional ways is an active research problem. 
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1 Introduction 



Since the discovery of higli-temperature superconductivity in the layered cuprates in 1986 [Bed- 
norz and Miiller 1986], it has been estabhshed that the superconducting order parameter has the 
d^2_y2-wa,vc symmetry in a number of these materials [Schrieffer 1994, Scalapino 1995, van Har- 
lingen 1995, Tsuei and Kirtley 2000]. This is important, because the d^2_y2-wave symmetry of the 
order parameter suggests the possibility of an electronically mediated pairing mechanism. Perhaps, 
the simplest model used for modelling the low-energy electronic correlations of the layered cuprates 
is the two-dimensional (2D) Hubbard model. Within this context, the nature of the pairing correla- 
tions in the Hubbard model as well as the nature of its low-lying electronic excitations has received 
considerable attention. 

In 1987 Anderson suggested that the 2D Hubbard model is relevant to the cuprates [Anderson 
1987]. However, even today questions remain about this model. In this article, what has been 
learned about the physical properties of the 2D and the 2-leg Hubbard models from the numerical 
studies will be reviewed. The emphasis will be placed on the d^2_y2 pairing correlations seen in 
these models and their microscopic origin. The implications of these calculations for the d^2_y2- 
wave pairing in the high- Tc cuprates is of current interest. In particular, one is interested in knowing 
whether the 1^2,2 _y2 -wave superconducting order could exist in the ground state of the 2D Hubbard 
model, and, if it does, whether it would have sufficient strength to explain the superconducting 
transition temperatures as high as those seen in the cuprates. There is also much interest in finding 
ways of enhancing the strength of the d^2_y2-wai:Ve pairing correlations observed in the Hubbard 
model. Within the past ten to fifteen years, the determinantal Quantum Monte Carlo (QMC) and 
the density-matrix renormalization-group (DMRG) techniques have been used to address these 
questions. These numerical studies are the subject of this review article. 

For the 2D Hubbard model, the QMC data on the magnetic, charge and the single-particle 
excitations will be presented. In order to investigate the pairing correlations, the QMC data on 
the irreducible particle-particle interaction and the solutions of the particle-particle Bethe-Salpeter 
equation will be shown. Using the DMRG method, the equal-time pair-field correlation function 
has been calculated in the ground state of the 2-leg Hubbard ladder. These DMRG data along with 
the QMC results on the 2-leg ladder will be shown and compared with each other. In addition, 
these results will be compared with the numerical studies of the t-J model. 

Figure 1.1 shows a schematic drawing of the CUO2 plane of the cuprates consisting of the one- 
electron Cu(3dj;2_j^2) and 0{px,Py) orbitals that give rise to the band in which the d^2_j^2-wave 
superconducting pairs form. In an electronic model of the Cu02 plane, there would be an onsite 
Coulomb repulsion at the Cu(3dx2_y2) and the 0(j)x,Py) orbitals and one-electron hopping matrix 
elements between the neighbouring Cu-0 and the 0-0 orbitals as well as longer-range hoppings 
and Coulomb interactions. A simplified approach is to use the single-band Hubbard model on the 
2D square lattice which has an onsite Coulomb repulsion and one-electron hopping matrix elements 
between the nearest-neighbour sites, as illustrated in Fig. 1.2(a). To some extent, this is motivated 
by the notion that the nonperturbative effects due to the onsite Coulomb repulsion will dominate 
at low temperatures and low energies. Extensive numerical calculations have been carried out for 
studying the properties of the one-band 2D Hubbard model. Another model which has proven useful 
for studying dx2_y2-wave pairing is the 2-leg Hubbard ladder illustrated in Fig. 1.2(b). This model 
has an onsite Coulomb repulsion U and the intrachain and interchain hopping matrix elements t and 
t±, respectively. The 2-leg Hubbard ladder is an important model where the pairing correlations 
can be studied in detail in the ground state for systems with up to 32 rungs. Here, the numerical 
results on the 2D and the 2-leg Hubbard models will be used to discuss the nature of the dx2_y2-wsLve 
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Figure 1.1: Schematic drawing of the Cu02 lattice consisting of the one-electron Cu{3d^2_y2) and 
the 0{px,Py) orbitals. 

pairing found in the Cu02 layers of the high-Tc cuprates. 

The Hubbard model was introduced for describing the Coulomb correlation effects in the d- 
bands of the transition metals, which show both band-like and atomic-like behaviour [Hubbard 
1963]. The Hubbard Hamiltonian has a one-body kinetic term and an interaction term which 
represents the onsite Coulomb repulsion when two electrons occupy the same orbital. In spite of 
studies covering four decades, questions remain about the Hubbard model. In two dimensions, the 
Hubbard Hamiltonian with only nearest-neighbor hopping is given by 

H = -t ^ {cl^Cja + C^j^Cia) + U^niani-^- H^Hia, (1.1) 

(i,j),a i ia 

where the sum over i and j is done over the nearest-neighbour sites on a square lattice. The one- 
electron hopping matrix element is t, the onsite Coulomb repulsion is U and the chemical potential 
/i is used for controlling the electron occupation in the grand canonical ensemble. Here, cj^ (cio-) 
creates (annihilates) an electron of spin a at site i and riia- = cj^Cicr is the occupation number of 
electrons with spin a at site i. In addition to the 2D Hubbard model, the numerical results on the 
2-leg Hubbard ladder with anisotropic hopping will be reviewed. In this case, the Hamiltonian is 

H = -tJ2 (4,A,<TCi+l,A,a + 4+l,A,aCi,A,a) " t± ^(4i,^Cj,2,a + 42,aCi,l,a) 
i,\,a i,(7 

where t is the hopping matrix element parallel to the chains and t± is the inter-chain hopping, as 
illustrated in Fig. 1.2(b). The operators c\ create an electron of spin a at site i of the A'th leg, 

and rii^x.a = c\ ^^Ci^x^„ is the electron number operator with spin a at site (i. A). 

The QMC and the DMRG are powerful numerical methods, which have been developed for 
studying strongly-correlated systems such as the Hubbard model. These techniques allow for the 
possibility to determine whether superconductivity exists in the ground state of an interacting 
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Figure 1.2: (a) Sketch of the 2D Hubbard lattice with isotropic nearest-neighbour hopping t and 
the onsite Coluomb repulsion U. (b) Sketch of the 2-leg Hubbard ladder with onsite Coulomb 
repulsions U, and the intrachain and interchain hopping matrix elements t and t±, respectively. 



system without first having to construct an approximate theoretical framework for describing its 
elementary excitations. The determinantal QMC method was introduced in Ref. [Blankenbecler 
et al. 1981]. The QMC results reviewed here were obtained by using the algorithm described in 
Ref. [White et al. 1989b]. Reviews of this method can be found in Refs. [Scalapino 1993, Muramatsu 
1999]. The DMRG method was developed by White [White 1992 and 1993], and since then has been 
used to calculate the equal-time correlation functions in the ground state of interacting systems. 
Here, in Sections 3 through 7, numerical results on the 2D and the 2-leg Hubbard models, which 
were obtained by using these two algorithms, will be presented. In Appendix, the determinantal 
QMC and the DMRG techniques will be described briefly. 

In Section 2, an introduction to the spin-fluctuation mediated d^2_y2-wa,ve superconductivity 
will be given. In addition, here, the experimental evidence for d^2_y2-wave superconductivity in 
the cuprates and the theoretical studies of d^2_y2-w8ive pairing in the cuprates will be reviewed 
briefly. 

In order to understand the microscopic many-body processes causing d^2_y2-wave pairing cor- 
relations in the 2D Hubbard model, it is useful to first discuss the magnetic, density and the 
single-particle excitations. The review of the numerical studies on the 2D Hubbard model will 
begin in Section 3 by presenting QMC data on the magnetic properties. Here, QMC results on 
the magnetic susceptibility will be shown at and near half-filling. In addition, an RPA-like model 
for the magnetic susceptibility will be introduced, and the effective particle-hole vertex will be 
discussed. The NMR experiments which support the existence of short-range AF fluctuations in 
the high-Tc cuprates will also be reviewed briefly. 

In Section 4, the QMC results on the charge susceptibility will be shown and its dependence 
on the strength of the Coulomb repulsion will be discussed. Here, the possibility of any ^^4]s.f" 
charge-density-wave (CDW) fluctuations in the 2D Hubbard model will be investigated. In this 
section, numerical results on the optical conductivity will also be discussed. 

The nature of the single-particle excitations in the 2D Hubbard model has received much atten- 
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tion within the context of the cuprates. In Section 5, the single-particle density of states N(uj) and 
spectral weight A(p,uj) obtained from the maximum-entropy analytic continuation of the QMC 
data on the single-particle Green's function will be reviewed. At half-filling, these calculations find 
a Mott-Hubbard gap, lower and upper Hubbard bands and quasiparticle bands. These quasiparti- 
cle bands are similar to those found in a spin-density-wave (SDW) insulator. It will be seen that, 
upon doping, a redistribution of the spectral weight takes place forming a narrow metallic band 
at the Fermi level. This band has unusually flat dispersion near the (vr, 0) and (0, tt) points of the 
Brillouin zone leading to large amount of single-particle spectral weight available for scatterings in 
the particle-particle channel. These results on A{p,uj) will also be compared with the results of 
the angular-resolved photoemission spectroscopy (ARPES) measurements on the cuprates. 

After reviewing the numerical results on the magnetic, density and the single-particle properties, 
in Section 6 the d^2_y2-wave pairing correlations in the 2D Hubbard model will be discussed. 
The QMC simulations have found that when doped away from half-filling there are short-range 
d^2_y2-wave pairing correlations in the 2D Hubbard model, but no long-range superconducting 
order has been observed in the parameter regime where the simulations are carried out. The 
irreducible particle-particle interaction Tj has been also calculated with the QMC simulations and 
it gives information on the microscopic many-body processes causing the attraction in the d^2_y2- 
wave pairing channel. In Section 6, these QMC results will be reviewed. It will be seen that 
Tj is repulsive at q = (7r,7r) momentum transfers and that it is the short-range AF correlations 
which are responsible for this behavior. Comparisons with the various diagrammatic approaches 
will also be carried out, which indicate that in the intermediate-coupling regime the momentum 
and the Matsubara-frequency structure in Tj can be described by a properly-renormalized single 
spin-fluctuation exchange interaction. Using the QMC data on Tj and the single-particle Green's 
functions, the Bethe-Salpeter equation in the particle-particle channel has been solved. These 
calculations allow for a comparison of the strength of the pairing correlations in the various channels. 
Here, it will be seen that as the temperature is lowered, the fastest growing pairing correlations 
occur in the singlet d^2_y2-wave channel. 

While these calculations are not carried out at sufficiently low temperatures to determine 
whether the d^2_y2-waYe long-range superconducting order exists in the ground state of the doped 
2D Hubbard model, the results on Tj are useful in the following sense. Consider a phonon-mediated 
superconductor such as Pb, where the effective interaction mediating the pairing forms already at 
temperatures of order the characteristic phonon frequency u>d- Hence, in this case, at T ~ lod it 
would be possible to study the effective particle-particle interaction responsible for superconduc- 
tivity, even though the superconducting long-range order takes place at much lower temperatures. 
Similarly, the QMC simulations are carried out at temperatures of order or less than the charac- 
teristic energy scale J ~ At^/U of the AF correlations. Hence, it is possible to probe Fj in the 
temperature regime where short-range AF correlations have formed, and see what type of many- 
body processes are important in mediating the pairing and which pairing channels are favored. 

The 2-leg Hubbard ladder is a model where the d^2_y2-waYe superconducting correlations can 
be studied in detail. In Section 7, the DMRG and the QMC results on the 2-leg Hubbard ladder 
will be reviewed. The DMRG calculations have found enhanced power-law decaying dr^2_y2-waYe 
superconducting correlations in the ground state of the 2-lcg Hubbard ladder. At the same time, 
with the QMC simulations, it is possible to calculate the magnetic susceptibility, the single-particle 
spectral weight and the particle-particle vertex Tj at sufficiently low temperatures. The compar- 
isons of the DMRG and the QMC results show that it is the strong short-range AF fluctuations 
which mediate the d^2_y2-wavc pairing correlations in this system. Furthermore, in this case, it 
is possible to study the dependence of the strength of the superconducting correlations on the 
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model parameters such as t_\_/t and U /t. It will be seen that the superconducting correlations are 
strongest in the intermediate-coupling regime and when t±/t is such that there is large amount 
of single-particle spectral weight pinned near the (tt, 0) and (0, tt) points of the Brillouin zone. In 
addition, the QMC calculations find that Vj peaks at (vr, vr) momentum transfer. These features of 
^(p, u) and F/ create optimum conditions for d,j.2 _y2-wsvc superconducting correlations. In Section 
7.3, these numerical results on the 2D and the 2-leg Hubbard models will be compared with each 
other. This comparison suggests that the d^2_y2-wave pairing correlations seen in these models do 
not require a particularly sharp peak in Tj but rather simply weight at momentum transfers near 
(vr,vr). 

In Section 8, the QMC and the DMRG results reviewed above will be compared with the 
diagrammatic and the other numerical approaches to the Hubbard model as well as with the 
numerical results on the t-J model. First, in Section 8.1, the QMC results on the 2D Hubbard 
model will be compared with the results of the fluctuation exchange (FLEX) approach. Here, 
the purpose will be to make simple estimates for the maximum possible in the 2D Hubbard 
model, if a superconducting phase were to occur. In particular, the effects of the system being near 
an AF Mott-Hubbard insulator on the strength of the d^2_y2-waYe pairing will be explored. In 
Section 8.2, the results from the variational and the projector Monte Carlo studies of the Hubbard 
model will be reviewed briefly. Recently, the dynamical cluster approximation (DCA) and the one- 
loop renormalization-group (RG) technique employing a 2D Fermi surface were applied to study 
d^2_y2 pairing in the 2D Hubbard model. The results of these calculations will be discussed in 
Section 8.3. The origin of the normal state properties of the cuprates in the pseudogap regime 
remains an important problem in this field. In Section 8.4, the findings of various calculations in 
the low-doping regime of the Hubbard model will be compared with the pseudogap regime of the 
cuprates. 

The Hubbard model is closely related to the t-J model, for which various Monte Carlo, exact 
diagonalization and DMRG calculations have been carried out. There is much interest in the 
ground-state phase diagram and the nature of the density and the pairing correlations in this model. 
In Section 8.5, the numerical studies on phase separation, stripe correlations and the d^2_y2-wave 
pairing correlations in the t-J model will be briefly reviewed and compared with the results on the 
Hubbard model. Comparisons will also be made with the 2-leg t-J model. 

In Section 8.6, the implications of the numerical results on the Hubbard model for the nature 
of the d^2_y2-wa:ve pairing in the high-Tg cuprates will be discussed. In Section 9, the summary 
and the conclusions will be given. 

2 da;2-y2-wave superconductivity in the cuprates 
2.1 Spin-fluctuation mediated d^2_y2-wave superconductivity 

Since the development of the BCS theory, it has been of interest to see whether the effective in- 
teraction which is responsible for pairing could be mediated by excitations other than phonons. 
The superfluidity of ^He is an example where it is believed that the pairing is due to the ex- 
change of ferromagnetic spin fluctuations resulting in triplet p-wave superconductivity [Leggett 
1975]. Within the context of organic superconductors, the possibility of the pairing being mediated 
by spin fluctuations had been noted in Ref. [Emery 1986]. For the heavy fermion materials, the 
possibility of d-wave superconductivity due to the exchange of AF spin flTictTiations was proposed in 
Refs. [Scalapino et al. 1986, Miyake et al. 1986, Cyrot 1986]. For the high-Tc superconductors, the 
possibility of d^2_y2-Yfme superconductivity was first studied in Ref. [Bickers et al. 1987], where 
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the framework of the two-dimensional Hubbard model was used near an SDW phase. 

Here, a discussion of these ideas on how the exchange of the AF spin-fluctuations leads to a 
d^2_y2-wave gap function will be given, since this picture will be useful for the remainder of the 
article. In the paramagnetic phase of the 2D Hubbard model and within RPA, the single spin- 
fluctuation exchange interaction between opposite-spin particles at zero frequency transfer is given 
by [Berk and Schrieffer 1966, Doniach and Engelsberg 1966] 

nPiP)-t/+ /W;'''^ + /y';'ri''. . (2.) 

l-f^Xo(P-P) l-t^Xo(p' + p) 

where Xo(q) is the usual Lindhard function. The longitudinal and the transverse spin fluctuations 
contributing to V are illustrated in Fig. 2.1. Here, T^(p'|p) is the interaction for the scattering of 
a pair of opposite-spin electrons at states (p, — p) to (p',— p')- Near half-filling and for a tight- 
binding band structure ep = —2t{cospx+cospy)—iJ,, where is the chemical potential, the Lindhard 
function peaks at q ~ (tt, tt). Consequently, for a system with Stoner-enhanced AF correlations, the 
spin-fluctuation exchange interaction F(p'|p) is large and repulsive at p' — p ~ (7r,7r) momentum 
transfers. Now, consider the BCS gap equation 

p' ^ 

where p and p' are restricted to being on the Fermi surface of the doped system as illustrated 
in Fig. 2.2. When this equation is solved, it is found that the leading superconducting instability 
occurs for a gap function with the (i^2_j^2-wave symmetry, 

Ap = ^ (cos - COS Py). (2.3) 

In order to see why the dx2_y2-wsNe form is a solution, consider a quasiparticle pair occupying 

(p, — p) with p near (vr, 0), as shown in Fig. 2.2. This pair can scatter to (p', — p') where p' ~ (0, vr) 
through the interaction l^(p'|p), which is enhanced and positive for momentum transfers near (tt, tt), 
because Ap changes sign between p and p'. This is basically the reason for why the dx2_y2-ws:ve 
gap symmetry is favored by the AF spin fluctuations. Since F(p'|p) is always positive, there is no 
singlet solution with the usTial s-wavc symmetry. For simplicity, these arguments were given at zero 
frequency, where p and p' are restricted to being on the Fermi surface, but the same arguments 
hold when the frequency dependencies of the gap equation and of the effective interaction are taken 
into account. In Section 6, the QMC data on the irreducible particle-particle interaction Vj and 
the solutions of the Bethe-Salpeter equation will be reviewed for the 2D Hubbard model. There, 
it will be seen that F/ is similar to the single spin-fluctuation exchange interaction and the leading 
singlet pairing instability occurs in the dr^2_y2-wme channel in the parameter regime where the 
QMC simulations are carried out. 



2.2 Experimental evidence for drc2_y2-^^^a^^e superconductivity in the cuprates 

Early experimental evidence for (i^2_y2-wave pairing in the cuprates arose from the measurement of 
the NMR longitudinal nuclear relaxation rate T^^ for Cu(2) and 0(2,3) nuclei in YBa2Cu307. In 
particular, the anisotropy of Tj"^ for Cu(2) below provided a signature for cia.2_j^2-wave pairing 
in this compound [Takigawa et al. 1991, Martindale et al. 1992, Bulut and Scalapino 1992]. These 
NMR measurements provided evidence for both the phase and the nodes of the c?^2_j^2-wave gap 



8 





Figure 2.1: Feynman diagrams illustrating the effective particle-particle interaction y(p'|p) within 
the single spin-fluctuation exchange approximation. 




Figure 2.2: Sketch of the Fermi surface of the 2D tight-binding model doped near half-filling. The 
two-particle momentum states (p, — p) and (p', — p') are indicated by the filled circles. 

function. The temperature dependence of the Cu(2) transverse nuclear relaxation rate was 
also used for the identification of the gap symmetry [Bulut and Scalapino 1991, Itoh et al. 1992]. 
These NMR measurements were also supported by the ARPES experiments which extracted the 
magnitude of the gap function |Ap| on the Fermi surface [Shen et al. 1993]. The microwave cavity 
measurements of the superconducting penetration depth A(r) in single crystals of YBa2Cu307 
found a linear temperature variation which also gave support to a d^2_y2-wave gap [Hardy et al. 
1993]. The phase-coherence measurements on YBCO-Pb dc SQUIDS [Wollman et al. 1993] and 
the tricrystal experiments of Tsuei et al. [Tsuei et al. 1994] established that in YBa2Cu307 the 
order parameter has the d^2_y2-wave symmetry. Reviews of these results on the symmetry of the 
gap function are given in Refs. [Schrieffer 1994, Scalapino 1995, van Harlingen 1995, Tsuei and 
Kirtley 2000] . Today, it is established that in a number of the hole and the electron doped cuprates 
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the gap function has the d^-i _y2-^f!&ye symmetry [Tsuei and Kirtley 2000]. 
2.3 Theoretical studies of (i^2_j/2-wave pairing in the cuprates 

For the high-T^ cuprates, the possibility of d^2_y2-^sye superconductivity was first studied by 
Bickers et al. [Bickers et al. 1987] using the framework of the 2D Hubbard model. In particular, 
the single spin-fluctuation exchange interaction was used for studying the d^2_y2-wave pairing near 
an SDW instability within the random-phase approximation. In these calculations, the effect of 
the single-particle self-energy corrections due to the spin fluctuations was taken into account. 

The next level of calculations for the dx2_y2-waYe superconductivity were carried out using 
the fluctuation-exchange approximation (FLEX) within the 2D Hubbard model [Bickers et al. 
1989]. This method self-consistently treats the fluctuations in the magnetic, density and the pairing 
channels. Within this approach it was found that the AF spin fluctuations lead to a d,j.2_y2-wave 
superconducting phase which neighbors the SDW phase. In these calculations, superconducting 
transition temperatures as high as 0.025t were found [Bickers et al. 1989, Bickers and White 1991]. 
Since the hopping matrix element within a one-band description of the cuprates is estimated to be 
about 0.45 eV [Hybertson et al. 1990], this value of Tc corresponds to about 130 K. However, the 
need to carry out exact calculations was also noted. 

Another approach to d^2_y2-wa.ye pairing in the cuprates was from a phenomenological point 
of view [Moriya et al. 1990]. In this approach, a spin-fluctuation exchange interaction was used to 
calculate the superconducting transition temperatures. The spin susceptibility and the parameters 
used in the model were obtained from fitting the NMR and the electrical resistivity data on the 
normal state of the cuprates within the self-consistent renormalization theory. In these calculations, 
the single-particle self-energy corrections due to the spin fluctuations were also included, and Tc's 
as high as those in the cuprates were obtained. A review of these calculations can be found in 
Ref. [Moriya and Ueda 2000]. 

A similar phenomenological approach was taken in Ref. [Monthoux et al. 1991]. In this ap- 
proach, the spin susceptibility used in the effective interaction was also taken from fltting the NMR 
data, and Tc's comparable to those in the cuprates were obtained. In these calculations, the Tc was 
initially calculated without including the self-energy effects. Later, the self-energy corrections were 
included self-consistently [Monthoux and Pines 1992]. 

These spin-fluctuation theories for superconductivity differ from phonon-mediated supercon- 
ductivity in a fundamental way. Here, the effective attractive interaction which is responsible for 
pairing the electrons is generated by the electronic correlations rather than by an external system 
such as the lattice vibrations. Hence, it is an important question whether in a microscopic model 
the electronic correlations can indeed produce an effective attractive interaction which leads to 
superconductivity, and, if so, what would be the nature and the strength of such pairing corre- 
lations. The QMC and the DMRG methods were developed to address these types of questions. 
The numerical studies of the Hubbard model carried out with this purpose are the subject of this 
review article. 

3 Magnetic correlations in the 2D Hubbard model 

The parent compounds of the cuprates are AF insulators, and the AF phase extends up to ^^5% 
doping. For dopings beyond this, there is clear evidence from the NMR and the inelastic neutron 
scattering experiments that the system exhibits short-range AF spin fluctuations [Pennington and 
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Slichtcr 1990, Birgeneau 1990]. These are properties which support using the Hubbard framework 
for studying the electronic correlations of the cuprates. 

The QMC calculations have found that the ground state of the Hubbard model at half-filling 
has AF long-range order [Hirsch and Tang 1989, White et al. 1989b]. Away from half-filling, the 
QMC calculations find short range AF correlations [Hirsch 1985]. In this section, the QMC results 
on the magnetic susceptibility x of the 2D Hubbard model will be reviewed, and compared with 
an RPA-like approach for modelling it. In addition, the QMC results on the irreducible particle- 
hole vertex will be discussed. These will be useful for gaining a microscopic understanding of the 
magnetic correlations in the Hubbard model. 




Figure 3.1: (a) Static magnetic susceptibility x(q, 0) versus q. (b) Equal-time correlation function 
](m^(r)m^(0))] versus r along (1, 0). These results are for U = 8t and (n) = 1.0 on an 8 x 8 lattice. 



3.1 Antiferromagnetic long-range order at half-filling 

The transverse magnetic susceptibility is defined by 

^('1''^) ^Jl^Jo dre'^"^-e-'^-'{m-{i + i,0)m+{i,0)). (3.1) 

Here, m+(i,0) = c||Cjj, and m~{i,T) = e^'^m^ {i,0)e~^'^ , where m~(i,0) is the hermitian conju- 
gate of m'^{i,0). Ref. [Mahan 1981] can be consulted for the finite-temperature Green's function 
formalism. In the following, x will be plotted in units of t~^. Figure 3.1(a) shows x(q, ia;^ = 0) 
versus q along various cuts in the Brillouin zone as the temperature is lowered. Here, it is seen that 
as T is lowered, xi^.^ 0) develops a sharp peak at the antiferromagnetic wave vector q = (tt, tt). 
Figure 3.1(b) shows the equal-time correlation function 

|(m-(r)m^(0))|, (3.2) 
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where m^(r) = cJ^|Cr| — c|!j^Crj^, versus r as the temperature is lowered on an 8 x 8 lattice at half- 
filling. In this figure, it is seen that as T is lowered below 0.5t, the AF correlation length S^af 
reaches the size of the system and long-range AF order is established on the 8x8 lattice. For an 
infinite lattice, the long-range AF order would take place at T = 0. 




Figure 3.2: (a) Static magnetic susceptibility xi^i 0) versus q. (b) Equal-time correlation function 
|(m^(r)m^(0))| versus r along (1, 0). These results are for U = 8t and (n) = 0.87 on an 8 x 8 lattice. 



3.2 Short-range antiferromagnetic fluctuations away from half-filling 

Figure 3.2(a) shows xi^i 0) versus q for (n) = 0.87 and U = 8t on the 8x8 lattice. Here, as T 
decreases, the development of a broad peak centred at q = (tt, tt) is seen. Figure 3.2(b) shows 
the corresponding |(m^(r)m^(0))| versus r. At the lowest temperature of 0.33t, the AF correlation 
length is about 0.6 lattice spacing for this filling. Here, ^ was obtained by fitting the decay of 
|(m^(r)m^(0))| to e'l-'l/^. 

The Monte Carlo results on xi^ii^m) have been analytically continued to the real frequency 
axis using the Fade approximation. The results on the spin-fluctuation spectral weight Imx(q, a;) 
versus a; at q = (tt, tt) obtained in this way are shown in Fig. 3.3(a). While the Fade approximation 
can only resolve the crude features of the spectrum, the softening of the AF spin fluctuations at 
low temperatures is clearly seen. 

Another quantity which can be used for characterising the low-frequency magnetic correlations 

is 

Imx(q,w) 



F(q) = "^"'^^ . (3.3) 

This quantity is of interest because it is used to calculate the longitudinal nuclear relaxation rates 
of the planar Cu and O nuclear spins in the layered cuprates. Figure 3.3(b) shows i^(q) versus q 
for (n) = 0.87 and U = 8t. As T is lowered, -F(q) gets enhanced at (7r,7r), and more weakly at 
q ~ 0. 
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Figure 3.3: (a) Spin-fluctuation spectral weight Imx(q, w) versus w at q = (vr,7r). (b) F{q) versus 
q for the same temperatures as in (a). These results are for U = 8t and (n) = 0.87 on an 8 x 8 
lattice. 

3.3 RPA-like model for the magnetic susceptibility 

An RPA-like approach has been used for modelling x(qj^) in the intermediate coupling regime 
with the form [Bulut et al. 1990, Bulut 1990, Chen et al. 1991, Bulut et al. 1993] 

X{(i,^) = - r (3.4) 

1 - ^xo(q,w) 

where xoi^j^) is the Lindhard function for the band electrons 



xo(,q,'^j 



1 ^ /(gp+q) - /(gp) . . 



p 



Ep = —2t{cospx + cospy) — fi and /(gp) is the usual Fermi factor. The parameter U entering 
Eq. ( |3.4| ) represents the reduced Coulomb repulsion, which is assumed to be renormalized due to 
the electronic correlations. In this simple approach, U is taken to be independent of momentum 
and temperature. This turns out to be a reasonable approximation in the temperature range where 
the Monte Carlo simulations have been carried out. 

Figure 3.4(a) compares this RPA-like approximation for x(q) 0) with the Monte Carlo data for 
U = 4t, T = 0.25t and (n) = 0.87. A similar comparison is given in Fig. 3.4(b) for the T dependence 
of x(q = (^)'^))0). In these figures U = 2t has been used. This choice will be discussed further in 
the following subsection, where results on the irreducible particle-hole vertex will be reviewed. 
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Figure 3.4: (a) Comparison of the QMC and the RPA results on x(q, 0) versus q at T = 0.25t. 
(b) Comparison of the QMC and the RPA results on x(Q = (7r,7r),0) versus T/t. Here, the QMC 
data (filled circles) are for U = 4t and the RPA results (solid curves) were obtained using U = 2t. 
In both figures, the results are shown for (n) = 0.87. 



3.4 Irreducible particle-hole interaction 

In the RPA form discussed above, a reduced Coulomb repulsion U which is independent of q is 
used. In order to test this approximation, the effective particle-hole vertex defined by 

U (q, ium) = 7 ^4 r, (3.6) 

X(q, lUJm) X(q, iOJm) 

was calculated in Ref. [Bulut et al. 1995]. Here, x(q, iuim) is the piece of the magnetic susceptibility 
which does not include the reducible particle-hole vertex, and it is defined by 

X(q,za;^) = -— ^ G{t[> + (\,iun + Wm)G{^,iUn), (3.7) 

where G(p, iu;„) is the exact single-particle Green's function calculated with QMC. In Fig. 3.5(a), 
the QMC data on x(q, 0) versus q is compared with the Lindhard function Xo(ct)O) at T = 0.25t. 
Here, one sees that x(q; 0) is suppressed with respect to Xo(q5 0), and this is because of the single- 
particle self-energy effects. The effective irreducible vertex L'"(q, 0) obtained from the QMC data 
on x(cii 0) s-iicl x(ct5 0) is plotted in Fig. 3.5(b) as a function of q at various temperatures. In this 
figure it is seen that U (q, 0) shows little dependence on q and T. This is the reason for the good 
agreement between the QMC data and the RPA results obtained using a reduced U in Eq. (3.4). 
However, the values of f/(q, 0) are larger than U = 2t used in the RPA. For instance, at T = 0.25t 
one has U{c[ = (7r,7r),0) = 3.2t. This is because in the RPA form the Lindhard susceptibility xo is 
used rather than x, which has the effects of the single-particle self-energy corrections. At this point. 
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Figure 3.5: (a) Comparison of the QMC data on x(q, 0) versus q with the Lindhard function 
Xo(q, 0) at T = 0.25t. (b) Effective particle-hole vertex ?7(q, 0) versus q along (1, 1). These results 
were obtained for U = At and (n) = 0.87 on an 8 x 8 lattice. 

it should be noted that the momentum dependence of the effective particle-hole vertex within the 
context of the cuprates was also discussed in Ref. [Anderson 1997]. 

With the QMC simulations the irreducible particle-hole vertex Fj was also calculated [Bulut et 
al. 1995]. The calculation of F/ is similar to that of the irreducible particle-particle vertex F/, which 
will be discussed in Section 6. In this approach, F/ is obtained from the solution of the particle-hole 
t-matrix equation illustrated in Fig. 3.6. These calculations showed that the QMC data on F/ are in 
agreement with the results on the effective particle-hole vertex U (q, 0) discussed in this section. In 
particular, F/ with the center-of-mass momentum q = (vr, vr) and frequency ujm = was averaged 
over the incoming and outgoing momenta near the Fermi surface, and it was found that this average 
(F/) is 3.6t ± 0.7t for U = 4t, T = 0.25t and (n) = 0.87. This agrees with C7(q = (7r,7r),0) = 3.2t 
found above. Hence, the effective particle-hole vertex for center-of-mass momentum q = (tt, tt) and 
frequency ujm = is about 3.2t for a bare Coulomb repulsion of U = 4t. The difference reflects the 
renormalization of the bare Coulomb repulsion in the particle-hole channel due to the higher-order 
many-body processes such as the Kanamori type of particle-particle scatterings [Kanamori 1963] . 
However, in the RPA form of Eq. (3.4), U = 2t needs to be used in order to take into account the 
single-particle self-energy corrections which are neglected within RPA. 

These results on the effective particle-hole vertex will also be useful when the irreducible particle- 
particle interaction F/ is compared with the single spin-fluctuation exchange interaction Tj^ in 
Section 6.4. There, it will be seen that when the effective coupling gll between the electrons 
and the spin fluctuations is taken to be about 3.2t, the resulting Tjf agrees with the QMC data 
on F/. One expects gll to be closely related to F/ and, hence, the QMC data on the effective 
particle-particle and particle-hole vertices are consistent. 
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Figure 3.6: Feynman diagrams for the particle-hole t-matrix equation. 



3.5 Comparison with the NMR experiments on the cuprates 

The measurements of the longitudinal NMR rate T^^ have provided evidence for the existence of 
short-range AF fluctuations in the normal state of the high- Tc cuprates. These experiments are 
reviewed in Refs. [Pennington and Slichter 1990, Takigawa 1990]. The relaxation rate Tf^ is given 



where |^(q)p is the hyperfine form factor of the particular nuclear spin, and F{q) has been defined 
in Section 3.2. The hyperfine form factors used in calculating the various NMR rates for the 
cuprates were derived in Ref. [Mila and Rice 1989]. 

An additional experimental quantity which supports the existence of the short-range AF cor- 



The measurements of for ^^Cu(2) in the normal state of the YBa2Cu307, YBa2Cu408 and the 
La2-a;Srj;Cu04 systems have provided evidence for the presence of the short-range AF correlations 
in these compounds. [Pennington and Slichter 1991, Itoh et al. 1992, Imai et al. 1993]. 

Since the Monte Carlo calculations cannot be carried out at low temperatures, it is difficult to 
make direct comparisons with the and the T2 measurements on the cuprates using the Monte 
Carlo data. However, the RPA-like model described above in Section 3.3 has been used to calculate 
the NMR rates in the normal state. Even though it is a simple approximation, this approach has 
been useful for analysing the T^^ and the Tg"^ data on YBa2Cu307 [Bulut et al. 1990, Bulut and 
Scalapino 1991]. These calculations imply that the low-energy AF spin fluctuations estimated for 
the 2D Hubbard model have the right order of magnitude for describing the general features of 
the AF correlations in optimally doped YBa2Cu307. The NMR data on YBa2Cu307 were also fit 
phenomenologically [Millis et al. 1990] and by using self-consistent renormalization theory [Moriya 
et al. 1990]. 

4 Charge fluctuations in the 2D Hubbard model 

In this section, the QMC data on the charge susceptibility and the optical conductivity of the 2D 
Hubbard model will be reviewed. 



by 




(3.8) 



q 
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4.1 Chcirge susceptibility 



The momentum and frequency dependent charge susceptibiUty is defined by 



n(q,za;^) 




dre' 



-S(q,r) 



(4.1) 



where 



^(q,r) = 



(rp(q,r)pt(q,0)), 



(4.2) 



with pT(q) = Cp_,_qo.Cpcr and p(q, r) = e p(q)e . The equal-time density-density 

correlation function S'(q) = S'(q, r = 0) and the charge compressibility k = d{n)/dn given by 
n(q 0,0) arc also of interest. In the following 11 and S will be plotted in units of t~^. The 
nature of the density correlations in the 2D Hubbard model was studied with various numerical 
techniques. The QMC [Moreo and Scalapino 1991, Moreo et al. 1991] and the exact diagonaliza- 
tion [Dagotto et al. 1992b] calculations did not find any indication of phase separation in the 2D 
Hubbard model. The charge compressibility of the Hubbard model near half-filling was calculated 
by zero-temperature QMC simulations [Furukawa and Imada 1991 and 1992]. At fillings (n) = 0.5 
and 0.2, 5'(q) was calculated with the determinantal QMC for U = At and 8t on up to 16 x 16 
lattices [Chen et al. 1994]. The Green's function Monte Carlo (GFMC) method was also used for 
calculating S(q) for the 2D Hubbard model on lattices with up to 162 sites in the ground state 
[Becca et al. 2000]. These calculations were carried out for U = 4t and lOt at and near half-filling. 
In this section, the QMC results on n(q, w) and S(q) from Ref. [Bulut 1996] will be shown for 
various values of U/t near half-filling. 

Extensive numerical calculations have been also carried out for studying the density correlations 
in the ground state of the t-J model. These calculations have found phase-separated, striped, and 
uniform phases, and they will be briefly described in Section 8. The results on n(q, w) and ^(q) 
presented here were obtained at relatively high temperatures and it is difficult to reach conclusions 
about the ground state of the Hubbard model using them. Nevertheless, here, the QMC data 
on the density correlations will be compared with the results on the t-J model obtained with the 
high-temperature series expansions [Putikka et al. 1994]. This gives valuable information on the 
dynamics of the density fluctuations in the 2D Hubbard model 

The high-temperature series expansions have found that the structure in S'(q) can be described 
by the noninteracting spinless-fermion model. This is due to the fact that the double-occupancy 
of the sites is not allowed in the t-J model, and once this constraint is taken into account, the 
density correlations are similar to those of the noninteracting particles. Consequently, for a filling 
{n), S{q) has structure at wave vectors which are associated with the Fermi surface corresponding 
to a filling of 2(n). In other words, the structure in S'(q) follows the "4ki7" wave vectors rather 
than the "2ki?" wave vectors of the system. 

First, results at half-filling will be presented. Figure 4.1 shows QMC results on n(q, 0) versus 
q at various temperatures for U = 8t. Here, n(q, 0) is rather featureless in q, and it gets further 
suppressed as T is lowered. This is due to the Mott-Hubbard gap which exists at half-filling. Next, 
results for {n) = 0.87 are shown. Figure 4.2(a) shows the temperature evolution of n(q, 0) versus 
q for [/ = 8t. In this case, the q ~ component of n(q, 0) gets enhanced as T decreases. The T 
dependence of n(q 0, 0) seen in Fig. 4.2(a) is consistent with the results of the zero-temperature 
QMC calculations of the charge compressibility [Furukawa and Imada 1991 and 1992]. Figure 4.2(b) 
shows Imn(q, w) versus u; for q = (7r/4, 0), where it is seen that the density fluctuations soften 
as T is lowered. The spectral weight ImH was obtained by the Fade analytic continuation of the 
QMC data. 
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Figure 4.1: Momentum dependence of n(q, 0) at half-filling for U = 8t and various temperatures. 

The evolution of the q structure in n(q, 0) with U/t is shown in Fig. 4.3(a) for T = 0.5t and 
(n) = 0.87. For comparison, in Fig. 4.3(b) results are shown for no(q, iw^ = 0) of the U = 
system given by 



In this figure, the solid curve represents no(q, 0) for (n) = 0.87 and the dotted curve is for (n) = 
1.74, which corresponds to n(q, 0) of the noninteracting spinlcss-fcrmion system for (n) = 0.87. 
As U jt increases, it is seen that n(q, 0) develops features which are more similar to those seen for 
the noninteracting spinless-fermion system. In Fig. 4.3(a), it is seen that the peak in n(q, 0) shifts 
from q = (7r,7r) to q ~ as ?7/t increases. It is not possible to obtain this type of momentum 
dependence for 11 (q, 0) from the simple RPA form 



where ?7 is a constant representing the effective coupling in the density channel. 

Next, QMC results on S'(q) will be discussed. Figure 4.4 shows <S'(q) versus q at various 
temperatures for U = 8t. Here, it is seen that S'(q) exhibits small variation with T/t. Figure 4.5(a) 
shows the evolution of S'(q) with U/t for T = OM and (n) = 0.87. Here, it is seen that with 
U/t increasing, S'(q) gets suppressed. However, in contrast to the case of n(q, 0), there is no 
qualitative change in the q structure of 5(q) as U/t increases, and the peak remains at (vr, vr). It 
is useful to compare these QMC results with >S'o(q) of the noninteracting system. Figure 4.5(b) 
shows S'o(q) versus q at T = O.bt for (n) = 0.87 (solid curve) and (n) = 1.74 (dotted curve). 
Comparing Figs. 4.5(a) and (b), one observes that the QMC data on ^(q) is more similar to <S'o(q) 
for (n) = 0.87, and the q structure in S{q) does not follow that of the noninteracting spinless- 
fermion system. These results on ^(q) are consistent with those found by the GFMC technique for 
the ground state [Becca et al. 2000] , and by the QMC method at finite temperatures [Chen et al. 
1994]. 
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Figure 4.2: (a) Momentum dependence of n(q, 0). (b) Imn(q,a;) versus w at q = (7r/4, 0). These 
results are for (n) = 0.87, U = 8t and various temperatures. 
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Figure 4.3: (a) n(q, 0) versus q at (n) = 0.87 for various values of U. (b) no(q, 0) versus q for 
the U = system at (n) = 0.87 (solid curve) and (n) = 1.74 (dotted curve). These results are for 
T = OM. 
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Figure 4.4: Equal-time density correlation function S'(q) versus q at (n) = 0.87 for U = 8t and 
various temperatures. 
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Figure 4.5: (a) ^(q) versus q at (n) = 0.87 for various values values of U. (b) <S'o(q) versus q of 
the U = system for {n) = 0.87 (solid curve) and {n) = 1.74 (dotted curve). These results are for 
T = 0.51 
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Here, it has been seen that the results on S{q) for the 2D Hubbard model are quite different 
from those calculated for the 2D t-J model for J ~ 0.5t with the high-temperature series expansion 
[Putikka et al. 1994] or by the GFMC technique [Calandra et al. 1998]. In the 2D Hubbard model 
near half-filling, S'(q) does not exhibit obvious features which might be associated with the 'Mkp" 
wave vectors of the Fermi surface. However, H(q, 0) has structure, which might be associated 
with the "4kj7" wave vectors, when U/t is large. Hence, the dynamical density fluctuations in 2D 
Hubbard model with large U/t have features which might be related to the 'Mkp" wave vectors, 
while the static density correlations appear to follow the "2ki;'" wave vectors. However, results on 
H(q, Laj) at lower temperatures and on bigger lattices arc necessary before reaching conclusions. 

The nature of the density correlations in the Hubbard model will be discussed further in the 
remainder of the article. In Section 7.3, comparisons will be made with the 2-leg Hubbard ladder, 
in which case power-law decaying "4kF" CDW correlations have been found [Noack et al. 1996]. 
In addition, in Section 8 there will be further comparisons with the density correlations in the 2D 
t-J model. 



4.2 Optical conductivity 

In this section, the results of the numerical calculations on the optical conductivity ai{u)) of the 
2D Hubbard model will be discussed. The optical conductivity for the 2D Hubbard model was 
calculated with the exact diagonalization technique for a 4 x 4 lattice with U = lOt at half-filling 
and in the doped case [Dagotto et al. 1992c]. In these calculations, the insulating gap in ai(uj) at 
half-filling is clearly seen. Upon doping to {n) = 0.875, the amount of the spectral weight above 
the insulating gap is reduced. In this case, the additional features found in ai^uo) are a Drude peak 
at a; = and spectral weight induced at intermediate energies. The Drude weight of the Hubbard 
model was also calculated with the QMC simulations [Scalapino et al. 1992 and 1993]. A review 
of the exact diagonalization calculations for cJi(ti;) and of comparisons with the experimental data 
on the cuprates can be found in Ref. [Dagotto 1994]. 

The real part of the frequency-dependent q = conductivity is given by 

(4.5) 

where the current-current correlation function A^x is defined as 

Kx{iu^m)= dTe~'^-^mx{r)j^m (4.6) 

with 

= -iteaY^{c\^Ci+xs - c\+xs'^is)- (4.7) 

is 

The analytic continuation in Eq. (4.5) is carried out using the maximum entropy procedure. In the 
following, the hopping t, the electron charge e and the lattice constant a are set to unity. 

The solid curve in Fig. 4.6 shows ai{u)) versus uj at half-filling for U = 8t, T = 0.125t and an 
8x8 lattice. These are data from Ref. [Bulut et al. 1994a]. Here, in spite of the limited resolution 
of the analytic continuation procedure, the insulating gap in ai{Lo) is clearly seen. The dotted curve 
in this graph is the mean-field SDW result given by 



cri{uj) = Re 



lUJr, 
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Figure 4.6: Real part of the frequency dependent conductivity of the half-filled Hubbard model 
with U = 8t (solid curve). These results were obtained for an 8 x 8 lattice at T = 0.125t. The 
dotted curve is the mean-field SDW result with A = 2At. A finite broadening of F = 0.5t has been 
used in plotting the SDW result. 

which was calculated using A = 2At. Near the threshold u ~ 2A, the SDW coherence factor 
^[1 - {el - A2)/e2] in Eq. 4.8 goes to 1 for £p near zero, and hence the peak at the threshold does 
not vanish within SDW. A check of the maximum entropy analytic continuation is the /-sum rule 
which has the form 

poo ^ 

dioai{io) = -{-k:c), (4.9) 

where (kx) is the average kinetic energy per site associated with hopping in the x direction. The 
area under the solid curve in Fig. 4.6 is 0.74, where the separate QMC measurement of ^{—kx) gives 
0.77. The difference reflects the difficulty in analytically continuing the numerical data. When the 
system is doped, it is difficult to estimate the accuracy of the analytic continuation procedure for 
(Ji((j), especially in the limit a; — > 0. Hence, away from half-filling, the maximum entropy or the 
Fade techniques were not used to extract (Ti{uj) or the dc resistivity p{T) from the QMC data on 
the current-current correlation function Axx- 

Even though it has not been possible to reliably extract p{T) for the 2D Hubbard model, for 
2D t-J clusters cJi(i^) and p{T) were calculated at various fillings by using the finite-temperature 
Lanczos technique [Jaklic and Frclovsck 2000]. For J/t = 0.3, the results for (Ti(uj) and p{T) were 
compared with the experimental data on cuprates in various doping regimes. In some cases, the 
agreement obtained with the experimental data was at a quantitative level. In the intermediate 
doping regime, where 0.75 < (n) < 0.85, it was shown that (Ti{uj) is not consistent with the usTial 
Drude form but rather with the marginal Fermi liquid concept where the only uj scale is given 
by T [Varma et al. 1989]. The finite-temperature Lanczos calculations find also that p{T) is 
proportional to T in this regime but the slope of p{T) changes at T ~ J. The underdoped and 
the overdoped regimes were also studied with this technique. Jaklic and Prelovsck concluded that 
the main features of the unusual normal state properties of the cuprates could be attributed to the 
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large degeneracy of states and the frustration induced by doping the antiferromagnet. 

5 Single-particle properties of the 2D Hubbard model 

The angular-resolved photoemission spectroscopy (ARPES) measurements on the high-Tc cuprates 
probe the singlc-particic properties of these compounds and have found a number of interesting 
features. In the normal state of the optimally doped cuprates, the ARPES experiments have found 
that there are quasiparticle-like bands which cross the Fermi level leading to a large Fermi surface. 
An unusual feature is that these bands have extended flat dispersion near the (7r,0) and (0, tt) 
points in the Brillouin zone [Dessau et al. 1993, Gofron et al. 1993]. The ARPES experiments 
have also provided valuable information aboT.it the evolution of the singlc-particic spectral weight 
with doping in the underdoped regime as the insulating state is approached. Reviews of these 
experiments can be found in Refs. [Shen and Dessau 1995, Damascelli et al. 2001]. Below, the 
numerical studies of the single-particle spectral weight in the 2D Hubbard model will be reviewed. 
In particular, the origin of the correlated metallic band which develops at the Fermi level upon 
doping the AF Mott-Hubbard insulator will be investigated. It will be seen that the 2D Hubbard 
model provides an explanation for some of the features found in the ARPES data. 

The single-particle properties of the 2D Hubbard model were studied with various many-body 
techniques. The single-particle spectral weight was calculated phenomenologically by taking into 
account the scattering of the quasiparticles by the AF spin fluctuations [Kampf and Schrieffer 1990a 
and 1990b]. The evolution of the single-particle spectral weight with doping was also studied with 
the exact diagonalization technique [Dagotto et al. 1991]. There have been extensive numerical 
calculations of ^(p, u;) by analytically continuing the QMC data on the imaginary-time single- 
particle Green's function [White et al. 1989c, White 1991, Vekic and White 1993, Bulut et al. 
1994a, Bulut and Scalapino 1995, Moreo et al. 1995, Haas et al. 1995, Duffy and Moreo 1995, 
Preuss et al. 1995 and 1997, Grober et al. 2000] . In these studies, the max;imum-entropy technique 
has been the main algorithm for the analytic continuation procedure. This technique was also used 
for calculating A{p,uj) for the three-band Cu02 model [Dopf et al. 1992b] and the one-dimensional 
Hubbard model [Preuss et al. 1994, Zacher et al. 1998]. These studies have provided valuable 
information about the single-particle spectral weight in the Hubbard model. In the following, the 
maximum-entropy results on the 2D Hubbard model will be reviewed. Here, the emphasis will be 
on the narrow correlated band which forms at the Fermi level upon doping the insulator, and on 
the fiat bands which are observed near the (tt, 0) and (0, vr) points in the Brillouin zone. 

The single-particle spectral weight ^(p, a;) is given by 

A{tp,uj) = ——lu\G{p,i(jOn ^ oj + i5) (5.1) 

TT 

where 

G{p,iu;n)= dTe'^-^Gii>,T) (5.2) 
Jo 

and 

G(p,r) = -(rcp.(r)ct^(0)). (5.3) 

With the maximum-entropy technique, one uses the QMC data on G(p, r) to invert the integral 
equation 

/oo e~^'^ 
-oo J- "r 6 
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in order to obtain ^(p,a;). By applying the same procedure to 



G,,{t) = -{TcUMM), (5.5) 

the single-particle density of states N{io) is obtained. The maximum-entropy technique for analyt- 
ically continuing QMC data is described in Refs. [Silver et al. 1990, Jarrell and Gubernatis 1996]. 
The statistical errors for Gii{T) are smaller than for G(p,r), hence the results for N{uj) shown 
below have higher resolution. Here, first the results on the evolution of N{uj) with doping will be 
reviewed, and later ^(p, to) will be discussed. In the following, N[lo) and A(p, cj) will be plotted 
in units of t~^. The data which will be shown in this section are from Refs. [Bulut et al. 1994a, 
Bulut and Scalapino 1995]. 




Figure 5.1: Single-particle density of states N(uj) versus u. These results were obtained for U = 8t 
and (n) = 1.0 on a 12 x 12 lattice. 



5.1 Single-particle density of states 

Figure 5.1 shows the single-particle density of states N{u) versus u> at half-filling for U = 8t on 
a 12 X 12 lattice at T = 0.5 and 0.125t. Here, the chemical potential /x is located at = 0. In 
this figure, it is seen that at T = 0.5t, N(uj) consists of two peaks corresponding to the lower 
and the upper Hubbard bands, which are separated by the Mott-Hubbard pseudogap. At this 
temperature, the AF correlation length is still less than the system size, as it was seen in Section 
3.1. As r is lowered down to 0.125t, the AF correlation length reaches the system size, and in 
this case additional sharp peaks appear at the upper edge of the lower Hubbard band and at the 
lower edge of the upper Hubbard band. At T = 0.125i, the pseudogap has become a full gap in the 
single-particle spectrum with a magnitude of 2 A ^ 4M. The size of the gap 2 A was also calculated 
within an SDW approach where the single-particle self-energy corrections were included [Schrieffer 
et al. 1989]. Within this approach, for U = 8t, 2A is found to be about 4.8t, while for U = 4t, 2A 
is about 2t. 
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The sharp peaks which are located below and above the Mott-Hubbard gap appear when the 
system has long-range AF order. These narrow bands have a bandwidth of about It, corresponding 
to « 2 J where J = At'^/U, and they exhibit SDW-like dispersion, which will be discussed in the 
next section. It is known that in the half-filled 2D t-J model the quasiparticle bandwidth is about 
2 J [Liu and ManoTisakis 1992]. The exact diagonalization calculations for the 2D t-J model also 
find SDW-like quasiparticle excitations with a bandwidth of about 2. 2 J at half- filling [Dagotto 
1994]. In the t-J model at half-filling, one doped hole propagates by flipping the spins in the AF 
background around it. For this reason, the bandwidth for the hole motion is determined by the 
magnetic exchange J rather than the hopping matrix clement t. The QMC data implies that the 
hole propagation is accompanied by similar many-body processes in the half-filled 2D Hubbard 
model. 




Figure 5.2: Single-particle density of states N{lo) versus uj at fillings (a) (n) = 0.94, (b) 0.87 and 
(c) 0.70. These results are for T = OM and U = 8t. 

Next, results on the evolution of N{oo) with doping are shown. Figure 5.2 shows A^(ti;) versus 
u at fillings (n) = 0.94, 0.87 and 0.70 for U = 8t and T = 0.5t. Upon hole doping, the chemical 
potential moves rapidly from a; = to the top of the lower Hubbard band. This is accompanied 
with a gradual transfer of spectral weight from above the Mott-Hubbard gap to the Fermi level. The 
transferred spectral weight goes into forming a narrow metallic band at the Fermi level. These data 
are in agreement with the exact diagonalization calculations [Dagotto et al. 1991]. In the overdoped 
regime at (n) = 0.70, the Fermi level lies below the metallic band. The temperature variation of 
N{u;) is shown in Fig. 5.3 for (n) = 0.87 and U = 8t. Here, the build up of a narrow band with 
a width between 4 J and 5 J is clearly seen as T decreases. The comparison of Figs. 5.3 and 5.1 
indicates that this correlated metallic band forms out of the SDW-like bands of the half-filled case. 

At this point, it should be noted that the general features of N{uj) seen in Fig. 5.3 are similar 
to what is found for the infinite-dimensional Hubbard model away from half- filling [Jarrell 1992]. 
N{uj) for the infinite-dimensional Hubbard model is calculated by first mapping the model to an 
Anderson impurity problem, which is then solved exactly. In this case, a narrow peak at the Fermi 
level is also found in addition to the lower and upper Hubbard bands. 

It is useful to compare these results on N(lo) for U = 8t with the results for U = At and 12t, 
which are shown in Fig. 5.4. For U = 12t, the narrow metallic band is further separated from the 
lower and the upper Hubbard bands, and the pseudogap is bigger. On the other hand, for U = At, 
the pseudogap is not observed at these temperatures, even though there are hump-like structures 
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Figure 5.3: Development of the single-particle density of states N{uj) as the temperature is lowered 
for (n) = 0.87 and U = 8t. 




Figure 5.4: N{io) versus u for (a) U = At and (b) U = 12t at (n) = 0.87 and T = OM. 

which might be attributed to the lower and the upper Hubbard bands. In addition, the peak at 
the Fermi level is broader. Hence, as U increases from 4t to 12t, the distribution of the weight 
in N{uj) changes considerably. However, the density of states at the Fermi level N{iJ,) varies by a 
small amount. 

It is also useful to discuss how N(iJ,) varies with doping at fixed U and T. Figure 5.5 shows 
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N{p,) versus (n) for U = 8t at T = 0.5t (filled circles) and l.Ot (open circles). It is clearly seen 
that at these temperatures A^(/i) has a peak at finite doping, and a depression in N(fx) exists near 
half-filling. However, it is not possible to extract the T ^ behavior. 
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Figure 5.5: Filling dependence of the density of states at the Fermi level Ar(/x) for U = 8t at 
temperatures T = OM (filled circles) and l.Ot (open circles). 

5.2 Single-particle spectral weight 

In this section, results on the single-particle spectral weight ^(p, to) will be discussed. It is useful 
to first study the effects of varying U on A{p,u}). Figure 5.6 shows A{p,aj) versus lo at various 
temperatures for U = 8t and 4t. These results were obtained for the p = (7r/2,7r/2) point on an 
8x8 lattice. For both values of U, a heavily-damped quasiparticle-like peak develops near the 
Fermi level as the temperature decreases. However, for U = 8t this peak is broader and it has 
reduced weight. In addition, in this case, the peak location varies with the temperature. Even 
though the maximum-entropy technique has an intrinsic broadening which is difficult to estimate, 
these results suggest that the damping of the quasiparticlcs is stronger for larger U/t. 

In the previous section, it was seen that for U = 8t a narrow metallic band develops at the Fermi 
level upon doping. In order to gain insight into the origin of this band, in Fig. 5.7, ^(p, lo) versus 
LU is plotted for p taken along various cuts in the Brillouin zone as indicated in the insets. These 
data were obtained on a 12 x 12 lattice for U = 8t, {n) = 0.87 and T= 0.5t. In this figure, spectral 
weight representing the upper and the lower Hubbard bands are seen in addition to a quasiparticle 
band which crosses the Fermi level. The quasiparticle peak is especially of interest. When p is on 
the Fermi surface, this peak is centred at w = 0. As p moves outside of the Fermi surface, the 
quasiparticle peak is clearly resolved. However, when p is below the Fermi surface, for instance at 
p = (0,0), the quasiparticle peak is obscured by the lower Hubbard band, even though the exact 
diagonalization calculations [Dagotto et al. 1992a] find a quasiparticle peak at this momentum in 
addition to the lower Hubbard band. This is because the maximum-entropy algorithm used here 
has poor resolution when u is away from the Fermi level. Indeed, similar calculations with higher 
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Figure 5.6: Comparison of ^(p = (7r/2, 7r/2), w) versus uj for (a) U = 8t and (b) U = At. Here 
results arc shown at various temperatures for (n) = 0.87, and the frequency axis is shifted such 
that the Fermi level occurs at a; = 0. 

resolution were able to observe the quasiparticle peak below the Fermi surface [Preuss et al. 1995]. 
An unresolved issue is whether the dispersing band near the Fermi level corresponds to a true 
quasiparticle band. However, the exact determination of the energy and the temperature scaling 
of the spectral weight cannot be carried out because the maximum-entropy technique has finite 
resolution and, in addition, bigger lattices and lower temperatures are required. 

An important feature of the quasiparticle band is that it is unusually flat near the (tt, 0) and 
(0, tt) points of the Brillouin zone. This causes large amount of spectral weight to be near the 
Fermi level. Even though the weight in the quasiparticle peak near the Fermi level decreases as 
U increases from At to 8t, the density of states at the Fermi level stays nearly the same. This is 
because of the flat bands which are pinned near the Fermi level. The presence of the flat bands is 
especially important because they generate phase space for the scattering of the quasiparticles in 
the d.j.2_y2-wave pairing channel. These features found in the 2D Hubbard model are similar to the 
ARPES data on the cuprates [Dessau et al. 1993, Gofron et al. 1993]. 

In addition to the QMC data, in Ref. [Dagotto et al. 1994], the relation with the fiat bands 
seen in the ARPES data were noted by using the dispersion of one hole in the half- filled t-J model 
calculated on a 16 x 16 lattice. It should also be noted that Beenen and Edwards have found good 
agreement with the main features of the QMC data using a two-pole approximation in calculating 
the single-particle Green's function [Beenen and Edwards 1995]. Dorneich et al. have extended 
this type of calculations for the Hubbard model for the case of (n) = 1 [Dorneich et al. 2000]. 
Finally, it should be noted that a number of transport measurements on the cuprates have been 
interpreted in terms of the Fermi level being close to a van Hove singularity [Tsuei et al. 1992, 
Newns et al. 1994]. Here, it is seen that in the 2D Hubbard model there can be extended flat bands 
near the Fermi level and they are produced by the many-body effects rather than being due to the 
one-electron band structure. 
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Figure 5.7: Single-particle spectral weight A{p.uS) versus lo for various momentum cuts through 
the Brillouin zone for U = 8t and (n) = 0.87. Here, the frequency axis is shifted such that the 
Fermi level occurs at a; = 0. 



These calculations for A{p, u) were also carried out at half-filling, and the results were fitted 
with the SDW form [Bulut et al. 1994a], 

A{p, u) = ul6{cj - Ep) + vlS{uj + Ep) (5.6) 

where = i(l + jp/Ep), = ^(1 - 'Jp/Ep), 7p = -2t{cospx + cos{py)), and Ep = ^Jll + A2, 
For large A/t, 

' " (5.7) 



= y 7i + A2 ~ A + J(cOSPa: + COSPyY. 

Similar calculations with higher resolution [Preuss et al. 1995] were carried out, which showed that, 
in fact, the spectrum is better described by ~ A + ^{cospx + cospyY, hence, the quasiparticle 
band at half-filling has a width of 2 J in agreement with the results on N(uj) at half-filling and with 
the calculations on the t-J model with one hole. 

Based on the results discussed above, schematic dispersion relations were constructed for the 
2D Hubbard model, which are shown in Fig. 5.8 for (n) = 0.87 and (n) = 1.0. In Fig. 5.8(a), the 
dashed horizontal line represents the chemical potential, and the solid curve denotes the narrow 
metallic band located at the Fermi level. The shaded regions represent the lower and the upper 
Hubbard bands. The width of the quasiparticle band is of order 4J. Figure 5.8(b) illustrates 
the distribution of the spectral weight at half-filling. Here, the solid and the dotted lines denote 
the SDW-like quasiparticle bands. The picture which emerges from these calculations is that the 
narrow metallic band which lies at the Fermi level for (n) = 0.87 forms out of the SDW-like bands 
which exist in the half-filled case. 

The QMC results seen in Fig. 5.7 were obtained at T = 0.5t for (n) = 0.87, where the AF 
correlation length ( is less than one lattice spacing. The QMC simulations were also carried out 
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at r = 0.25t for (n) = 0.93, where ^ is about 1.2 lattice spacing [Preuss et al. 1997]. In this case, 
spectral weight has been observed at wave vectors which are displaced by Q = (vr, vr) with respect 
to the narrow metallic band at the Fermi level. These results are consistent with the calculations 
where " shadow bands" due to the scattering of the quasiparticles by the AF spin fluctuations were 
found [Kampf and Schrieffer 1990b]. This type of shadow structures were also observed in the 
photoeniission experiments [Aebi et al. 1995]. Detailed single-particle spectra in the low-doping 
regime of the 2D Hubbard model were obtained recently [Grober et al. 2000]. 




(0.0) (tt.tt) (tt.O) (0.0) (0.0) (tt.tt) (tt.O) (0.0) 

P P 

(a) (b) 

Figure 5.8: Schematic drawing of Ep versus p at J7 = 8i for (a) (n) = 0.87 and (b) (n) = 1.0. 
Here, the shaded areas represent the lower and the upper Hubbard bands, and the solid and the 
dotted curves indicate the quasiparticle bands. In (a), the solid points indicate the position of the 
quasiparticle peak obtained from Fig. 5.7 when the quasiparticle peak is not obscured by the lower 
Hubbard band, and the horizontal dashed line denotes the chemical potential ji. 

In this section, it was seen that the single-particle properties in the 2D Hubbard model are 
determined by the AF spin fluctuations and the Coulomb correlations. An especially important 
feature of the single-particle spectral weight is the flat bands near p = (7r,0) and (0,7r). They 
create large amount of phase space for the scattering of the quasiparticles in the d^2_y2-w&ye 
pairing channel. In Section 7, it will be seen that similar flat bands exist in the 2-leg Hubbard 
ladder, and they play a key role in determining the strength of the d^2_j^2-like superconducting 
correlations in this system. In Section 8.1, these QMC results on the single-particle spectral weight 
will be compared with those obtained from the FLEX approximation, and the implications for the 
dj;2_j^2-wave pairing in the 2D Hubbard model will be discussed. 
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6 Pairing correlations in the 2D Hubbard model 



An important question for the high-Tc cuprates is whether the 2D Hubbard model exhibits super- 
conductivity, and if does, what is the nature of the pairing. Various views were expressed on this 
subject, and some of them are as follows. In one view, the effective interaction between the particles 
is dominated by the exchange of an AF spin fluctuation, and this leads to dj.2_j^2-wave pairing [Bick- 
ers et al. 1987]. In an alternative view, the magnetic correlations are responsible for the pairing, 
but it is because of the "bag effect" which leads to pairing with the extended s or the d^2_y2-waYe 
symmetry [Schricffcr et al. 1988 and 1989]. In this approach, the effective pairing interaction is 
qualitatively different than the single spin-fluctuation exchange interaction. Another view is that in 
the ground state of the doped 2D Hubbard model there are no single-particle excitations carrying 
both charge and spin and the one-layer Hubbard model does not exhibit superconductivity [Ander- 
son 1987, Anderson et al. 1987, Anderson and Zou 1988]. In this view, the pairing is mediated by 
interlayer hopping. In order to differentiate among these theories, QMC simulations were carried 
out for the 2D Hubbard model. In this section, these QMC calculations will be reviewed. 

In Section 6.1, the QMC results on the pair-field susceptibilities will be reviewed. These calcu- 
lations have found that there is an effective attractive interaction in the singlet dx2_y2-wave channel 
[White et al. 1989a]. In addition, it is found that there is an attractive interaction in the extended 
s-wave channel. However, at the temperatures where the simulations are carried out, the pairing 
correlations are short range and do not show scaling with the lattice size [Moreo and Scalapino 
1991, Moreo 1992]. QMC simulations were also carried out for the pair-field susceptibilities in the 
three-band Cu02 model and similar results were found [Scalcttar et al. 1991, Dopf et al. 1992a]. 

Using QMC simulations, the irreducible particle-particle interaction Tj was calculated [Bulut 
et al. 1993, 1994b, 1995] and these results will be reviewed in Section 6.2. In these calculations, 
it is found that the momentum and the Matsubara-frequency structure in Tj follows closely that 
of the magnetic susceptibility, which means that the AF spin fluctuations dominate the effective 
interaction in the parameter regime where the simulations are carried out. In Section 6.3, the 
Bethe-Salpeter equation in the particle-particle channel will be solved using the QMC data on Fj 
and the single-particle Green's function G. The solution of the Bethe-Salpeter equation makes it 
possible to determine the strength of the various pairing channels quantitatively. Finite-size scaling 
of these results shows that as T is lowered, the fastest growing pairing instability is in the singlet 
dj.2_j^2-wave channel. Here, it is also found that as U/t increases from 4 to 8, the strength of 
the d,j.2_y2 pairing correlations grows. Later in Section 8.1, these calculations will be compared 
with the results of the FLEX approximation. In Section 6.4, the results on Fj will be compared 
with the single spin-fluctuation exchange approximation. In Section 6.5, comparisons will be made 
with the perturbation theory results for Fj, which are third order in U. These comparisons are 
useful for gaining insight into the effects of the various subgroups of many-body scattering processes 
contributing to Fj. 

6.1 Pair-field susceptibilities 

In order to study the pairing correlations in the 2D Hubbard model, the pair-fleld susceptibilities 
defined by 



where Aq is a pair creation operator with a symmetry, were calculated. The QMC calculations of 
the pair-field susceptibilities for the Hubbard model were first carried out in Ref. [Hirsch 1985]. In 




(6.1) 



31 




Figure 6.1: Feynman diagrams illustrating the pair-field susceptibilities (a) Pa and (b) Pq. Here, 
r represents the reducible particle-particle interaction. 

these calculations, it was found that there is an attractive effective interaction between two electrons 
with antiparallel spins when they are separated by one lattice spacing. In this reference, the pair- 
field susceptibilities in the singlet s and extended s-wave (s*) channels and the triplet channels were 
calculated, however the singlet d^2_y2-wave channel was not considered. Later in Ref. [White et al. 
1989a], various pair-field susceptibilities including the one with the singlet dx2_y2-wa,ve symmetry 
were calculated with QMC. In this section, results from Ref. [White et al. 1989a] will be shown for 
Pa which were calculated using the following pair-field operators in the dj.2_y2, s, extended s (s*) 
and p-wave channels, 

Ai =li2i-^y4+iAv (6-2) 
At =4^4, (6.3) 

Al* =\E4^eAv (6.4) 
At = -(ct_^ ^ct, -ct .ct,). (6.5) 

Here i sums over the four neighbours of site i. 

In order to see whether a given pairing channel with the symmetry a is attractive. Pa was 
compared with Pa, where Pa is the component of Pa which does not include the reducible particle- 
particle interaction T as illustrated in Fig. 6.1. Figure 6.2 compares Pa and Pa as a function of 
T/t for U = At and (n) = 0.87 on an 8 x 8 lattice. Here, it is seen that Pa is enhanced with respect 
to Pa by the particle-particle reducible vertex for the d,j.2_y2-waYe and weakly for the extended 
s-wave symmetries. On the other hand, in the s- and j;-wave channels the pair-field susceptibility 
is suppressed with respect to Pa- This means that the effective particle-particle interaction is 
attractive in the d^2_y2 and s* channels, and it is repulsive in the s and p-wave channels. However, 
at the lowest temperatures where the QMC simulations can be carried out, the d^2_y2-waYQ pairing 
correlations are only short range: P^ does not grow as the system size increases. Furthermore, 
while Pd is enhanced with respect to Pd, it is suppressed with respect to P^^, the d.-^2_y2-wawe 
pair-field susceptibility of the C/ = system, as seen in Fig. 6.2(a). So, at these temperatures, 
the dx2_y2-wave pair-field susceptibility gets suppressed when the Coulomb repulsion is turned on, 
which is not encouraging for d^2_y2-wave superconductivity in the 2D Hubbard model. 
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Figure 6.2: Pair-field susceptibifities Pa and Pa versus T/t for the d^2_y2, s*, s and p-wave pairing 
channels. These results are for U = At and (n) = 0.87 on an 8 x 8 lattice. In (a), Pq denotes the 
d^2_y2-waYe pair-field susceptibility of the U = system. From [White et al. 1989a]. 
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Figure 6.3: Feynman diagrams illustrating the correlation function A{j/ ,k'\p,k) in terms of the 
single-particle Green's functions G and the reducible particle-particle interaction T. 



6.2 Irreducible particle-particle interaction 

In this section, the QMC data on the irreducible particle-particle vertex Tj from Rcf. [Bulut et 
al. 1993, 1994b, 1995] will be reviewed. Using QMC simulations it is possible to calculate the 
two-particle Green's function 

A(a;4,X3|x2,a;i) = -(Tcf (a;4)q(a;3)c|(x2)c|(xi)), (6.6) 

where c}j{xi) with Xi = (xj,ri) creates an electron with spin a at site Xj and imaginary time Tj. By 
Fourier transforming on both the space and the imaginary-time variables one obtains 

A{p',k'\p,k) = -6pp:Skk'G^ip)G^{k) + ^Sk',p+k-p'G^{p')G^{k')r{p',k'\p,k)G^{p)G^{k), (6.7) 

where p = (p,za;„), Ga{p) is the single-particle Green's function, and r{p',k'\p,k) is the reducible 
particle-particle vertex. Hence, using the QMC data on A and G, T can be calculated. This 
equation is illustrated in Fig. 6.3 in terms of the Feynman diagrams. Here, the particle-particle 
interaction will be studied in the zero center-of-mass momentum and energy channel, and hence k 
and k' will be set to —p and —p', respectively. At the next stage, the irreducible particle-particle 
vertex Tj is obtained from the Monte Carlo results on F and G by solving the particle-particle 
i-matrix equation, 

Ti{p'\p) = Tip'lp) + ^Y1 riip'\k)G^{k)G^i-k)mp), (6-8) 

k 

which is illustrated in Fig. 6.4. In Eq. (6.8), T{p'\p) is used as a short notation for T{p', —p'\p, —p). 
This procedure for calculating F/ is essentially the opposite of the usual diagrammatic approach 
in which Fj is used to solve for F. In solving the t-matrix equation, an upper frequency cut-off of 
order the bandwidth is used. Finally, the singlet component of the irreducible vertex is obtained 
from 

^\vj(p'\p) + Vi{-p'\p) . 



(6.9) 



In the following, the QMC data on F/^ and F^ will be plotted in units of t as a function of the 
momentum transfer q = p' — p and of the Matsubara- frequency transfer cOm = ^n' — ^n- Here, p 
will be kept fixed at (vr, 0) and cOn at vrT. 

It is useful to start the discussion by first showing results on the singlet reducible vertex Fg. 
In Fig. 6.5, Fs(q, iw^ = 0) versus q is shown for U = At and 8t. As the temperature is lowered, 
the q ~ (tt, vr) component of F^ increases, and the q ~ component is suppressed. Here, it is seen 
that Tg becomes quite large especially for U = 8t. In Fig. 6.5(b), the q = (0,0) point is not shown 
because of large error bars for this point due to the way F is calculated from A. It is also not 
possible to show results at lower temperatures since the error bars grow rapidly. For instance, at 
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Figure 6.4: Feynman diagrams illustrating the i-matrix equation which relates the irreducible 
particle-particle interaction F/ to the reducible particle-particle interaction F. 
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Figure 6.5: Reducible particle-particle scattering vertex in the singlet channel Ts{q,,i(jJm = 0) 
versus q for (a) U = 4t and (b) U = 8t. Here, q = p' — p and p is kept fixed at (tt, 0). These 
results were obtained at (n) = 0.87 on an 8 x 8 lattice. For U = 8t, Fs(q, 0) at q = (0,0) is not 
shown because of large error bars in the data in this case. 

T = 0.33t and U = 8t, Ts{ci = (vr, tt), 0) was calculated to be 60i ± 20t after long simulation times, 
while Fsiq, = (7r/4, 7r/4), 0) was lOt ± lOt. 

Next, in Fig. 6.6(a) and (b), the momentum and the Matsubara-frequcncy dependence of the 
irreducible particle-particle vertex r/s(q, iuJm) is shown for U = At. At q = (vr, vr) momentum trans- 
fer, F/s(q, 0) reaches values larger than twice the bandwidth. It is useful to compare Figs. 6.6(a) 
and 6.5(a) in order to understand the effect of the repeated particle-particle scatterings on F/^. 
At these temperatures, the effect of these scatterings is basically to suppress the momentum and 
frequency-independent background in F/^. For instance, at T = 0.25t the difference in the magni- 
tude of F7s(q, 0) between q = (tt, tt) and (7r/4, 7r/4) is about lOi, which is the same as in Fs(q, 0). 
This is similar to the suppression of the screened Coulomb repulsion, which varies more slowly in 
frequency compared to the phonon propagator, in the usual phonon-mediated pairing. As a d^2_y2- 
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Figure 6.6: (a) Momentum and (b) the Matsubara-frequency dependence of the irreducible particle- 
particle scattering vertex in the singlet channel r75(q, iLUm)- In (a), ^is versus q = p' — p is plotted 

where p is kept fixed at (vr, 0). In (b), F/^ versus ujm is plotted for q = (vr, vr). Here, uJm = ujn' ~ 
and ujn is kept fixed at ttT. These results are for U = At and (n) = 0.87 on an 8 x 8 lattice. 

wave superconducting instability is approached, the expected behaviour for an infinite system is 
that Fs(q = (tt, 7r),0) +oo while Fs(q = (0,0), 0) — — oo. These QMC data show that such 
resonant scattering in the d,j.2 _y2-weive channel is not taking place at these temperatures. 

It is desirable to know how Fjs varies as U/t increases. However, it has not been possible to 
obtain Tjg for U = 8t from the t-matrix equation because of the larger error bars for this case. 
Nevertheless, in Fig. 6.5(b) it is seen that F5(q, 0) for U = 8t and T = OM exhibits large variation 
of order 20t between points (7r/4, 7r/4) and (vr, vr). If at these temperatures the eS'ect of the t-matrix 
scattering is only to suppress the q and tOm independent background in F/^, then this means that 
Tjs grows considerably as U/t increases from 4 to 8. 

Next, the temperature evolution of F7s(q, iujm) for ?7 = 4i is compared with that of the magnetic 
susceptibility xiHi'^-^m)- Figure 6.7(a) shows Monte Carlo data on x{^.^0) versus q at the same 
temperatures as in Fig. 6.6. The Matsubara frequency dependence of x(q = (vr, vr), iwm) is shown 
in Fig. 6.7(b). Comparing Figs. 6.6 and 6.7, one sees that the temperature evolution of F/s(q, i(jJm) 
closely follows that of x(q, ^t^^m)- Both of these quantities peak at q = (7r,7r), and as Um increases 
Fjs goes to the bare U value while x decays to zero. The relation between Tjg and the spin 
fluctuations will be studied in more detail in Section 6.4, where the QMC results on F/^ will be 
compared with the single spin-fluctuation exchange interaction. 

At this point, it should be noted that the q dependence of the effective pairing interaction for 
the cuprates was also studied within the spin-bag approach where a dip at q ~ is found in addition 
to the peak at q ~ (tt, vr) [Kampf and Schrieffer 1990a]. The dip is due to the cross-exchange of 
two AF spin fluctuations and it is responsible for the effective attraction between two spin bags. 
At r = 0.25i, which is the lowest temperature where the QMC calculation of Tjg can be carried 
out, a dip in Fjs at q ~ is not found. 
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Figure 6.7: (a) Momentum and (b) the Matsubara- frequency dependence of the magnetic suscep- 
tibiUty x{c\^,iijjm) at various temperatures. In (a), x{<\^Wm) versus q is plotted for = 0. In (b), 
x(q = (tt, tt), iwm) versus LOm is plotted. These results are for U = and (ra) = 0.87 on an 8 x 8 
lattice. 
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Figure 6.8: Real-space pattern of r/s(R) for U = 4t, T = 0.25t and (n) = 0.87. Here, it is seen 
that for R = 0, r/s(R) is strongly repulsive while being attractive when the singlet electron pair 
is separated by one lattice spacing. 
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In order to gain insight into the real-space structure of the effective particle-particle interaction, 
it is useful to consider the Fourier transform 

r,,(R) = ^Y1 e'^P~P'>^TUp',iu;n'\l>,ioJn) (6.10) 
p,p' 

for the lowest Matsubara frequencies uJn = ^n' = ^^T. Figure 6.8 shows the Monte Carlo data 
on r/s(R) as a function of R for T = 0.25t. At R = 0, F/^ is strongly repulsive, as expected, 
but for a singlet electron pair separated by one lattice spacing, F/^ is attractive. As the pair 
separation increases further, F/^ oscillates in sign and its magnitude decreases rapidly, reflecting 
the short-range nature of the interaction. Pairing correlations with the proper space-time structure 
can avoid the large onsite repulsion while taking advantage of the near-neighbour attraction. Thus, 
the interaction F/^ is attractive in the d^2_y2-wave channel. 

In Fig. 6.8, it is seen that the effective attractive interaction is —0.5t at the nearest-neighbour 
site. Here, one might argue that the long-range Coulomb repulsion between the electrons, which is 
not taken into account in the Hubbard model, could overcome this attraction. However, one would 
expect the long-range Coulomb repulsion to have weak uj^ dependence compared to that seen in 
Fig. 6.6(b), and after the t-matrix scatterings it should have weak influence on the strength of the 
pairing in the d^2_y2-wave channel. 

6.3 Bethe-Salpeter equation in the particle-particle channel 

In this section, the particle-particle Bethe-Salpeter equation will be solved using the QMC data on 
the irreducible interaction Tj and the one-electron Green's function. This way one can determine 
the magnitude of the eigenvalues and the (p, iun) structure of the leading pair-field eigenfunctions. 
This approach is useful for studying the leading scattering channels in the t-matrix quantitatively. 
For instance, one can determine how close the system is to a Kosterlitz-Thouless superconducting 
transition, and compare the strength of the pairing in various channels. The particle-particle 
Bethe-Salpeter equation is 

KMp) = ^iiP\p')G^ip')Gi{-p')Mp'), (6.11) 

P' 

where (paip) is the pair-field eigenfunction and is the corresponding eigenvalue. The Feynman 
diagram representing Eq. ( |6.11| ) is shown in Fig. 6.9. When the leading eigenvalue reaches one, 
the superconducting transition takes place, and the Bethe-Salpeter equation for the corresponding 
eigenfunction becomes equivalent to the superconducting gap equation. 




Figure 6.9: Feynman diagram representing the Bethe-Salpeter equation. 

Here, data from Refs. [Bulut et al. 1993] will be shown. Both the triplet and the singlet solutions 
of the Bethe-Salpeter equation will be considered. The momentum and the Matsubara-frequency 
structures of the leading eigenfunctions will be presented. In addition to the singlet (i^2_y2-wave 
eigenfunction (/)rf(p, zw), solutions which are odd in iVn and have s and p-wave symmetries are 
found. The odd-frequency eigenfunctions with s-wave and p-wave symmetries correspond to triplet 
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and singlet solutions, respectively. The finite-size scaling of the eigenvalues indicates that as T is 
lowered the singlet d^2_^2-wave eigenvalue grows fastest, and at low T the dominant singlet pairing 
channel has the dri.2 _y2-wme symmetry. In these calculations, the even-frequency extended s-wave 
channel was not found to be one of the leading pairing channels. 

The case of the odd-frequency pairing channel is interesting. The possibility of pairing in the 
triplet odd-frequency s-wave channel was first studied within the context of '^He in Ref. [Berezinskii 
1974]. The singlet odd-frequency p- wave channel was studied initially by using an effective interac- 
tion which is mediated by phonons [Balatsky and Abrahams 1992]. The possibility that an effective 
attraction in the singlet odd-frequency p-wave channel could be generated by the spin fluctuations 
was noted in Ref. [Abrahams et al. 1993] at about the same time the QMC calculations reviewed 
here were carried out. The general properties of the odd- frequency superconductors were discussed 
in Ref. [Abrahams et al. 1995]. 

In addition, here the effects of increasing U /t on the leading pairing channels will be discussed. 
As seen in the previous section, it was not possible to calculate Tig for U = 8t. However, if F is used 
rather than Tj in Eq. (6.11), then the corresponding eigenvalues are given by Aa/(1 — Aq,), where 
Aq, are the eigenvalues for the irreducible vertex. Hence, through this indirect way it is possible to 
study the leading Bethe-Salpeter eigenvalues for U = 8t. It will be shown that as U/t increases from 
4 to 8, the leading eigenvalues grow, since the effective particle-particle interaction gets stronger, 
but the momentum and the frequency structure of their eigenfunctions does not exhibit qualitative 
changes. 

In general, the Bethe-Salpeter equation can have both singlet and triplet solutions corre- 
sponding to a pair-wave function that has overall even or odd parity when p = (p, goes 
to (— p, —itjJn)- Here, the pair wave functions are characterized by its symmetry in momentum and 
spin space. The usual singlet s and d^2_y2-waYe states are even in frequency and even in momen- 
tum, (/)(p, —ibJn) = 4'{v-,'^'^n) and 0(— p,itj„) = 0(p, while the usual triplet or {py) state is 
even in frequency and odd when goes to —px- There are also odd- frequency pair- wave functions. 
In this case, one can have an odd-frequency s-wave triplet for which (/)(p, —itOn) = — '/'(p^^'^n) and 
4>{—p, iojn) = (j){p, iujn), or an odd-frequency px (or py)-wave singlet with (j){p, —ioJn) = — 0(p, iijJn) 
and (j){—p,iujn) = —(j){p,iuJn)- In the regime of the Hubbard model that is being studied, the 
s-wave triplet, and the p and (i^2„y2-wave singlet solutions are dominant. 

The temperature evolution of the four largest eigenvalues for U = At and (n) = 0.87 are given in 
Table 1. The momentum and the frequency dependence of the corresponding pair- wave functions 
are shown in Figs. 6.10 and 6.11 for T = 0.5t. At this temperature, an s-wave triplet state has 
the largest eigenvalue A^ 0.23. As seen in Figs. 6.10(a) and 6.11(a) (solid circles), the pair-wave 
function (/)s(p, iu;„) of the s-wave triplet state is even in p and odd in a;„. The open circles in 
Figs. 6.10(a) and 6.11(a) represent the pair-wave function (/)s/(p, «u;„) which has the second largest 
eigenvalue. This is also an odd-frequency s-wave triplet state. The dx2_y2-wme singlet state shown 
as the solid circles in Figs. 6.10(b) and 6.11(b) has the third largest eigenvalue. The fourth largest 
eigenvalue corresponds to a state which is odd in both ojn and p, having py (or px) symmetry [open 
circles in Figs. 6.10(b) and 6.11(b)], hence it is also a singlet. 

As T is lowered from l.Ot down to 0.25t, the largest eigenvalue Xg stays nearly the same, while 
\(i increases by about a factor of seven. This can be understood in terms of the temperature 
dependence of Ti{p'\p) which enters the Bethe-Salpeter equation. The pair- wave functions which 
are smooth in p but odd in oJn make optimum use of the (a;„, a;„') frequency structure of the repulsive 
Ti{p'\p) for pairing. However, as the temperature is lowered and Ti[p'\p) for p' — p = (vr, vr) grows, 
the dx2_y2 and p-wave solutions can make better use of the momentum structure in F/, and their 
eigenvalues get enhanced. 
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Figure 6.10: Momentum dependence of the leading Bethe-Salpeter eigenfunctions ^a(p,ia;„). 
These results are for a;„ = ttT, U = At, T = 0.5t and (n) = 0.87 on an 8 x 8 lattice. In (a), 
two odd-frcqucncy s-wave eigenfunctions are plotted, and in (b) the d^2_y2-w8ive (filled circles) and 
Py-wecve (open circles) are plotted. 



Q. 




Figure 6.11: Matsubara frequency dependence of the leading Bethe-Salpeter eigenfunctions 
4>a{p,iijJn)- These results are for U = At, T = OM and {n) = 0.87 on an 8 x 8 lattice. In 

(a), (p{p = (tt, 0),iu;„) is shown for the odd-6<J„ s and s' channels. In (b), (p{p,iujn) is shown for the 
singlet cia.2_2^2-wave channel at p = (tt, 0) and for the singlet py-wave channel at p = (7r/2, 7r/2). 
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T/t 


A. 


As' 


Xp 


Xd 


1.0 


0.261 


0.037 


0.023 


0.022 


0.50 


0.228 


0.090 


0.066 


0.076 


0.33 


0.251 


0.104 


0.095 


0.130 


0.25 


0.264 ± 0.007 


0.124 ±0.006 


0.129 ± 0.007 


0.182 ± 0.006 



Table 1: Temperature dependence of the Bethe-Salpeter eigenvalues for U = At and (n) = 0.87 on 
an 8 X 8 lattice. The error bars represent the uncertainty due to the Monte Carlo sampling. The 

error bars were not calculated at all temperatures, since it requires considerably more computer 
time. When the error bars are not indicated, they are estimated to be less than 10%. 



LxL 


A. 


Xp 


Ad 


4x4 
8x8 


0.296 ±0.010 
0.264 ± 0.007 


0.204 ±0.014 
0.129 ±0.007 


0.184 ±0.016 
0.182 ±0.006 



Table 2: Finite-size dependence of the Bethe-Salpeter eigenvalues for T = 0.25i, U = At and 
(n) = 0.87. 




Figure 6.12: Comparison of the dj.2_j^2-wave eigenfunction ^^{p,mT) with the usual (i^2_j^2-wave 
form Ap = (Ao/2)(cospj, ± cospy) for U = At, T = 0.25t and (n) = 0.87. Here, both (f)d{p,iTrT) 
and Ap have been normalised to 1 at p = (tt, 0). 



Table 2 shows the finite-size effects on the leading eigenvalues at T = 0.25i. Here, it is seen that 
the finite-size effects are especially large for the p-wave channel, and Ap decreases as the system 
size grows from 4 x 4 to 8 x 8. The finite size effects for the dx2_y2-wsive case are small. Hence, 
these results show that as T is lowered, grows fastest and at low T the dominant singlet pairing 
channel has the d^i-yi-wave symmetry. 

It is useful to compare the momentum dependence of the da;2_j^2-wave eigenvalue cpdip, iujn) at 
iOn = ttT with the usual (ij.2_j^2-wave gap form Ap = { Ao / 2) {cos Px — cos py), since this is often used 
in modelling the superconducting state of the cuprates. Figure 6.12 shows that the QMC data on 
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Figure 6.13: Temperature dependence of the d^2_y2-wave eigenvalue for U = At and 8t at 
{n) = 0.87. 

(j)ij_{p,iTTT) taken at T = 0.25t follow the usual d^2_y2-wsye form closely. 

Finally, the filled circles in Fig. 6.13 show the growth of \d for U = 4t as the temperature is 
lowered, while the open circles represent A,^ for U = 8i. The results on A^ for U = 8t were obtained 
as described above. Here, it is seen that A^ increases with U/t. 

At the lowest temperature where A^ can be calculated, the system is far from a Kosterlitz- 
Thouless superconducting transition which would be signalled hy ^ 1. Hence, while these data 
imply that at the temperatures where the QMC simulations are carried out, the singlet dj.2_y2-waYe 
pairing correlations are becoming dominant, it is not known whether A^ — 1 at lower temperatures. 
Below in Section 8.1, these results on A^ will be compared with the results of the FLEX calculations. 

The QMC calculation of A^^ for the doped case cannot be carried out at any lower temperatures. 
However, for half- filling it is possible to calculate the various pairing eigenvalues at low T and 
compare them with the eigenvalues of the Bethe-Salpeter equation in the AF particle-hole channel, 
Xa- The leading magnetic eigenvalue Ai occurs in an even- frequency s-wave channel with center- 
of-mass momentum Q = (vr,7r). In Fig. 6.14, Ai is compared with the various pairing eigenvalues 
for U = At and 8t at half-filling. As expected, at low temperatures Ai reaches 1 asymptotically, 
signalling the phase transition to the AF order state on the 8x8 lattice. Here, it is also seen that as 
the AF correlations develop at half-filling, A^^ becomes the leading pairing eigenvalue while always 
staying smaller than the magnetic eigenvalue Ai. 

6.4 Comparison with the spin-fluctuation exchange approximation 

Various spin- fluctuation exchange theories have been used for studying dj.2_y2-wave pairing for the 
cuprates [Bickers et al. 1987 and 1989, Moriya et al. 1990, Monthoux et al. 1991]. In this context, 
it is of interest to sec to what extent the Monte Carlo results for the irreducible vertex can be 
modelled by a single spin-fluctuation exchange interaction. For this purpose, here we compare the 
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Figure 6.14: Eigenvalues of the particle-hole and the particle-particle Bethe-Salpeter equations 
versus T for (a) U = At and (b) C/ = 8t for an 8 x 8 half-filled lattice. The solid points (•) are for 
the leading eigenvalue Ai of the Bethe-Salpeter equation in the AF particle-hole channel with the 
center-of-mass momentum Q = (vTjVr). The open symbols denote the even-frequency d^2_j^2-wave 
(o), and the odd-frequency p (A) and s-wave (□) eigenvalues of the particle-particle Bethe-Salpeter 
equation. 



Monte Carlo data with the approximate form 

rf^(q, iu;^) = U + ^g^U^ x(q, iu^m). (6.12) 

This form is motivated by the single spin-fluctuation exchange interaction [Berk and Schrieffer 1966, 
Doniach and Engelsberg 1966], which basically has this form with g = 1 near the antiferromagnetic 
instability. The Feynman diagrams illustrating the single spin-fluctuation exchange interaction 
were shown in Fig. 2.1. The factor of 3/2 arises from the two transverse and one longitudinal spin 
fluctuations. In calculating T^^ with Eq. ( |6.12| ), we will use Monte Carlo results for xi^i^i^m) and 
also set g = 0.8. The corresponding value of 3.2t for the effective coupling gU is consistent with 
the results of the Monte Carlo calculations of the effective irreducible vertex in the particle-hole 



channel, C/(q, 0), which were discussed in Section 3.4. Formally, Eq. ( 6.12 ) is analogous to the 
effective interaction in the electron-phonon superconductor 

V{q^,iuJm) =^ + Xl \9cixf Dx{q^,iuJm), (6.13) 
A 

where Dx{(i,iuJm) is the dressed phonon propagator and |(7qAp is the renormalised electron-phonon 
coupling. 

Figure 6.15(a) shows the single spin-fluctuation interaction in the singlet channel Tj^ (q, zcum) 
versus q at various temperatures. These results compare well with Tjs{(l,iujm = 0) seen in 
Fig. 6.6(a). Likewise, the comparison of Fig. 6.15(b) with Fig. 6.6(b) shows that the frequency 
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Figure 6.15: Single-spin-fluctuation exchange approximation for the irreducible particle-particle 
scattering vertex in the singlet channel Tis{(l,iuJm)- In (a), Tjg versus q is plotted. In (b), Tjs 
versus uj^ is plotted for q = (tTjTt). These results have been obtained for U = At and (n) = 0.87 
using the Monte Carlo data on x(q;^'^m)- 



dependence of Tj^ is in agreement with the Monte Carlo data. Considering the simplicity of 
Eq. ( |6.12 ), this agreement with the Monte Carlo data is quite good. These comparisons suggest 
that a properly renormalized single-spin-fluctuation exchange interaction is capable of reproducing 
the basic features of the effective particle-particle interaction in the weak-to-intermediate coupling 
Hubbard model at temperatures greater or of order J/2. 



6.5 Comparison with the perturbation theory 

In this section, the perturbation theory results for Tjg through third order in U from Ref. [Bulut 
et al. 1995] will be shown. These provide insight into the various subprocesses contributing to 
Tig. Figure 6.16 shows the diagrams contributing to Tjs up to third order in U. The dashed 
line in this figure represents the bare Coulomb repulsion. The diagrams (a) and (h) represent the 
first two terms in an RPA series corresponding to the exchange of a longitudinal spin fluctuation. 
Similarly, the low-order contributions arising from the exchange of a transverse spin fluctuation are 
represented by the diagrams (b) and (e). The diagrams (c) and (d) can be considered as corrections 
to diagram (b), where the bare particle- hole irreducible vertex is renormalized through Kanamori 
type of particle-particle scatterings [Kanamori 1963]. The diagrams (f) and (g) represent vertex 
corrections to the bare interaction (a). 

These diagrams have been evaluated on an 8 x 8 lattice with U = 4t, and the results of the 
various contributions are shown as a function of q for = in Fig. 6.17(a). The frequency 
dependence for q = (vr, vr) is shown similarly in Fig. 6.17(b). Here, the results on F/^ are plotted in 
the same way as in the previous section. In these figures, the filled circles represent the diagrams 
(a) and (h) in Fig. 6.16, and the open circles are for the diagrams (b) and (e). The triangles and the 
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Figure 6.16: Feynman diagrams contributing to the irreducible particle-particle interaction through 
third order in U. 



squares in Figs. 6.17(a) and (b) represent the contributions of the (c)-(d) and (f)-(g), respectively. 
The dominant contributions arise from the leading contributions to the spin-fluctuation exchange 
processes (a), (b), (e) and (h), and they give contributions which peak at q = (tt, tt). For lo^ = 0, 
both the Kanamori renormalization graphs and the vertex corrections act to reduce the strength 
of the q = (tt, tt) interaction by about one-third. In Fig. 6.18, the results obtained by summing 
the diagrams shown in Fig. 6.16 are compared with the Monte Carlo results for T = 0.25t. At this 
temperature, the low-order graphs in Fig. 6.16 are not adequate to represent the large momentum 
behaviour of the effective interaction. 

The importance of the vertex corrections to the single spin-fluctuation exchange interaction was 
pointed out in Ref. [Schrieffer 1995]. In particular, it was noted that for large values of the ratio 
x(q = (tt, tt), 0)/x(q 0,0), the vertex corrections should suppress the spin-fluctuation exchange 
interaction at q = (7r,7r) momentum transfer. The diagrams (f) and (g) in Fig. 6.16 are the lowest 
order vertex corrections, and, indeed, at this level they act to suppress the q = (tt, tt) component 
of F/s as seen in Fig. 6.17. However, the QMC data indicates that the suppression of the peak at 
q = (tt, tt) in Tjg is partial. This is because the ratio x(q = (tt, tt), 0)/x(q 0, 0) is about 5 rather 
than being of order 100 at the lowest temperature Fjs was calculated with QMC. Hence, a large 
AF correlation length is not necessary, and simply weight in Tjs at large momentum transfers is 
sufficient to yield a sizeable attractive interaction in the d.j.2_y2-wave channel. 

The results reviewed in this section constitute what has been learned about d^2_y2-wave pairing 
in the 2D Hubbard model from the determinantal QMC simulations. These simulations are not 
carried out at lower temperatures because of the sign problem, and, hence, it is not possible to know 
whether superconducting long-range order develops at lower T [Loh et al. 1990]. In Fig. 6.19(a), 
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Figure 6.17: Various diagrams contributing to the irreducible particle-particle interaction 
Tis{q_,iuOm) through third order in U. In (a) the momentum dependence is shown for Um = 0, 
and in (b) the frequency dependence is shown for q = (tt, tt) . Here, the filled circles represent the 
contribution of the bare U and the longitudinal spin-fluctuation exchange (diagrams of type (a) 
and (h) in Fig. 6.16), the open circles represent the contribution of the transverse spin fluctuations 
(Fig. 6.16(b) and (e)), the open squares show the ordinary vertex corrections (Fig. 6.16(f) and (g)), 
and the open triangles show the Kanamori type of vertex corrections (Fig. 6.16(c) and (d)). 

the average sign of the fcrmion determinants, (sign), which is defined in the Appendix, is plotted 
as a function of T at (n) = 0.87 for U = 8t and 4t. As the value of (sign) decreases below 1, the 
statistical error in the QMC data grows rapidly requiring exponentially long simulation times. In 
Fig. 6.19(b), the filling dependence of (sign) is shown for U = 8t at T = OM and 0.33t. These 
figures show the boundary of the parameter regime of the Hubbard model which cannot be probed 
because of the sign problem. However, the DMRG calculations [White 1992] are carried out at zero 
temperature, and they have provided valuable information about this regime in the 2-leg Hubbard 
ladder. The DMRG studies of the 2-leg Hubbard ladder will be reviewed in the next section. 

7 2-leg Hubbard ladder 

The 2-leg Hubbard ladder has been studied using various many-body techniques such as the DMRG 
[Noack et al. 1994, 1995, 1996 and 1997], the QMC [Dahm and Scalapino 1997], the exact diag- 
onalization [Yamaji and Shimoi 1994], and the weak-coupling renormalization group [Balents and 
Fisher 1996]. The results of the DMRG and the QMC calculations will be reviewed here. 

At half-filling, the ground state of this system does not have AF long-range order, rather it has a 
spin gap. Near half-filling, this model exhibits short-range AF correlations and power-law decaying 
dx2-y2-\ike superconducting and "4kir" CDW correlations [Noack et al. 1994 and 1996]. The 
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(a) (b) 

Figure 6.18: Comparison of the perturbation theory results on Tjs{c[,iu>m) through third order 
in U with the QMC data. In (a) the momentum dependence is shown for = 0, and in (b) the 
frequency dependence is shown for q = (tt, tt). 




Figure 6.19: (a) Temperature and (b) the filling dependence of (sign). In (a), the results are 
shown for U = 8t and At at (n) = 0.87, and in (b) for U = 8t. 
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DMRG calculations have found that in the ground state of this system the pair-field correlations 
for the interacting system can get enhanced with respect to those of the noninteracting (U = 0) 
system for a range of the model parameters. This is the first time ever in an exact ground state 
calculation for a bulk system that by turning on an onsite Coulomb repulsion the superconducting 
correlations get enhanced. For this reason, the 2-lcg Hubbard model is quite important. 

Another reason for studying this model is that here it is possible to understand the mechanism 
mediating the d^2_y2-\ike superconducting correlations by comparing the DMRG results with the 
QMC data obtained at relatively low temperatures. These comparisons indicate that it is the short- 
range AF fluctuations which mediate the pairing. Furthermore, the pairing is strongest when the 
model parameters are such that there is enhanced single-particle spectral weight near the (tt, 0) and 
the (0, tt) points of the Brillouin zone. For U = 8t and (ra) = 0.875, this occurs when t±/t ~ 1.5. In 
this case, the irreducible particle-particle vertex peaks at momentum transfers near (tt, tt) creating 
optimum conditions for d^2_y2 pairing. 

It is also interesting to study the 2-lcg Hubbard ladder because the half-filled insulating state 
is spin gapped while in the 2D case there is long-range AF order. In the cuprates, on the other 
hand, the undoped system has long-range AF order, and a spin gapped phase lies between the 
superconducting and the insulating phases. The 2-leg Hubbard model is a system where the relation 
between the spin gap and the superconducting correlations as well as the density correlations can 
be studied exactly. 

In section 7.1 below, the DMRG results on the pair-field correlation function from Ref. [Noack 
et al. 1997] will be shown for various values of the model parameters. In order to understand 
these results better, in Section 7.2 the QMC results on the single-particle spectral weight from 
Ref. [Noack et al. 1997], the irreducible particle-particle interaction and the solution of the Bethe- 
Salpeter equation from Ref. [Dahm and Scalapino 1997] will be presented. In Section 7.3, the results 
on the 2-leg Hubbard ladder will be compared with those on the 2D case. Later, in Section 8.5.3, 
comparisons will be made with the superconducting correlations found in the 2-leg t-J ladder. 

7.1 DMRG results 

The DMRG calculations found that there are power-law decaying dj.2_y2--wave pair-field correlations 
in the ground state of the 2-leg Hubbard ladder. The d^2_y2-like internal structure of the pairs can 
be seen by considering the pair-creation amplitude 

{N2\{cl^cl^-cl^cl,^)\N,) (7.1) 

for adding a singlet pair on near- neighbor sites along and across the legs. Here, \Ni) is the ground 
state with four holes relative to the half-filled band and \N2) is the ground state with two holes on 
a 2 X 16 ladder. The results on the pair amplitude are shown in Fig. 7.1 for U/t = 8 and t±/t = 1.5. 
Note the d^2_y2-like change in the sign of this matrix element. 

Using the DMRG method it is possible to calculate the rung-rung correlation function 

D{i,j) = {A{i)A\j)) (7.2) 

in the ground state of the 2-leg Hubbard ladder with the open boundary conditions. Here, 

creates a singlet pair across the z'th rung, and c|^^ is the electron creation operator with spin 
a at the i'th site of the fc'th leg of the 2-leg ladder. The fact that the matrix elements shown 
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Figure 7.1: Schematic drawing of the pair- wave function showing the values of the off-diagonal 

matrix element (-^2|(Cr|c|!,j^ — c|!j^c|!,|)|A/'i) for creating a singlet pair between near-neighbor sites. 
Here, it is seen that the matrix element for this process is negative when the singlet-pair is created 
along a chain, while it is positive across a rung. This shows the dx2-y2-like nature of the pairing 
correlations in the 2-leg Hubbard ladder. 

in Fig. 7.1 arc finite means that this bare pair-creation operator, which is composed of the bare 
electron-creation operators c|^^, has finite overlap with the true pair-creation eigen-operator for 
this system. In order to minimize the effects of the boundaries, here D(i,j) is averaged over six 

pairs with I = \i — j\ fixed. This averaging starts with symmetrically placed values and 
then proceeds to shift these to the left and right of the center. By comparing results obtained on 
different size lattices, it is possible to control the finite size effects. In the following, D{i) calculated 
on the 2 X 32 lattice for i < 20 will be shown. In this case, the finite size effects are negligible. 

Figure 7.2 shows D{(.) versus £ for various values of t±/t with U = 8t and (n) = 0.875. The 
dashed and the dotted lines represent power-law decays of £~'^ and £~^, respectively. In this figure, 
it is seen that D{1) exhibits a power-law decay for t±/t < 1.6. For t±/t = 1.0, D{i) decays as 

When t±/t is increased from 1.0 to 1.4, the strength of D{£) gets enhanced and it decays 
more slowly. For t±/t = 1.6, D{£) is reduced and it decays faster. Below in Section 7.1, it will 
be seen that when t±/t > 1.4, the antibonding single-particle band becomes unoccupied, and the 
decrease in the strength of the pairing correlations is due to this effect. Using the data in Fig. 7.2 
for 1 < £ < 18, D[£) has been fitted to a form 



with a linear least-squares approximation. The resulting 9 values are plotted as a function of t±/t 
for various fillings in Fig. 7.3. Here, it is seen that the minimum in 9 versus t±/t depends on the 
filling, but it occurs for t±/t ~ 1.4 near half-filling. For t±/t > 1.4, the antibonding band is no 
longer occupied and the pairing correlations decrease rapidly. Thus the pairing correlations are 
enhanced near the point at which the antibonding band moves through the Fermi level [Noack et 
al. 1995 and 1997, Yamaji and Shimoi 1994]. 

Another measure of the strength of the pair-field correlations which can be used is the average 




(7.4) 
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Figure 7.2: Pair-field correlation function D{t} versus i for various values oit±/t with U = 8t and 
(n) = 0.875. 



CO 



1.5 



0.5 



0.5 




2x32 
U/t=8 



• {n)=0.875 
o {n)=0.9375 

J , , \ , L. 



1 1.5 



Figure 7.3: Exponent 9 versus t±/t iov U/t = 8 and various values of (n). 



of D{i)/D{l) for rung separations £ = 8 to 12: 



1 



(7.5) 



Figure 7.4 shows D versus t±/t for U = 8t at various fillings. This clearly shows how sensitively 
the pairing correlations depend on the value of t±/t. Next, the variation of D with U/t is shown 
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in Fig. 7.5. In this figure, the crosses represent the results for U = 0. Hence, the onsite Coulomb 
repulsion can significantly enhance the pairing correlations for a range of t±/t values. Here, one 
also observes that the enhancement of D is strongest for U/t between 3 and 8, which is in the 
intermediate coupling regime. In addition, as U/t increases the value oft±/t at which the peak in 
D occurs shifts towards smaller values. Hence, the strength of the pairing correlations is a sensitive 
function of t±/t, U/t and (n). 




0.5 1 1.5 2 



tx/t 

Figure 7.4: Averaged pair- field correlation function D versus t±/t foi U = 8t and various values 
of (n). 

These DMRG results have a special place in many-body physics. They represent the first case 
in an exact ground-state calculation for a bulk system where the pair-field correlation function 
gets enhanced by turning on an onsite Coulomb repulsion. These calculations were carried out 
on ladders with up to 32 rungs resulting in negligible finite-size effects. As seen in Section 6.1, in 
the 2D Hubbard model, it was found that by turning on the Coulomb repulsion the d^2_y2-w&\e 
pair-field susceptibility gets enhanced with respect to P^. However, was always found to be 
suppressed with respect to the pair-field susceptibility P^ of the C/ = system. This meant that, at 
these temperatures, the effective attractive interaction in the d.j.2_y2-wave channel is not sufficiently 
strong to overcome the suppression of Pd induced by the single-particle self-energy effects. Here, it 
is seen that for the 2-leg Hubbard ladder in the ground state, D can get enhanced over the U = 
result for a set of the model parameters. 

7.2 QMC results 

In order to understand these DMRG results better, in this section QMC data on the 2-leg Hubbard 
ladder will be shown. 

7.2.1 Single-particle spectral weight 

The single-particle properties of the 2-leg Hubbard ladder at half-filling were studied in Ref. [Endres 
et al. 1996]. Here, the evolution of the single-particle spectral weight ^(k,a;) with t±/t will be 
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Figure 7.5: Averaged pair-field correlation function D versus t^jt for (n) = 0.9375 and various 
values of IJ jt. 



shown for the doped case. These data, which were obtained for a 2 x 16 ladder with periodic 
boundary conditions along the chains, are from Rcf. [Noack ei al. 1997]. The comparison of the 
results on yl(k,a;) and D{£) will show that the pairing correlations are enhanced for t±/t such that 
the bottom of the antibonding band moves through the Fermi level. 

The results on A{]s.,uj) were obtained by the maximum-entropy analytic continuation of the 
Monte Carlo data. Figures 7.6 and 7.7 show data for U/t = 2 and 4, respectively. In these figures 
(n) = 0.94 and T = 0.125t, and the results are given for various values of tj^/t. Here, ^(k,a;) 
is shown for both k± = Q (bonding) and k± = tt (antibonding) bands as a contour plot in the 
uj-k plane where the intensity of the shading represents the magnitude of y4(k, cj) and k is the 
momentum along the chains. In the C/ = system, the quasiparticle dispersion consists of the 
bonding and the antibonding bands given by 

ek = -2icosA;±2i_L, (7.6) 

where k = {k,k±). The dotted curves in Fig. 7.6 represent for U = 0. Here, it is seen that 
the dispersion obtained for U/t = 2 closely follows that of the U = system. Note also that 
for t±/t = 1.8 the bottom of the antibonding band is located right at the Fermi level while for 
t±/t = 2.0 the antibonding band becomes unoccupied. In Fig. 7.5, it was seen that the pairing 
correlations for U/t = 2 are strongest when t±/t = 1.8. This comparison suggests that it is the 
variation in the single-particle spectral weight with t±/t which controls the dependence of the 
pairing correlations on t±/t. 

Figure 7.7 shows similar data on A(k,a;) for U/t = 4 and t±/t = 1.4, 1.6 and 1.8. In this case, 
there arc more differences between the QMC data and the U = results denoted by the dotted 
curves. Here, the antibonding band becomes unoccupied for a smaller value of t±/t = 1.8. One 
also notices that the top of the bonding and the bottom of the antibonding bands are flattened, 
especially for t±/t = 1.6, increasing the amount of the single-particle spectral weight near the Fermi 
level. This behaviour is similar to the build up of spectral weight near (tt, 0) in the 2D Hubbard 
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(a) 




(b) 




(c) 




Figure 7.6: Distribution of the single-particle spectral weight A{k, to) in the k and lo plane. The 
intensity of the shading indicates the amount of the spectral weight. These results are for U = 2t, 
{n) = 0.94, T = 0.125t and (a) t±/t = 1.6, (b) 1.8 and (c) 2.0. 



model as discussed in Section 5, and in the ARPES experiments on the high- Tc cuprates which are 
reviewed by [Shen and Dessau 1995]. 

In the next section, it will be seen that for the 2-leg Hubbard ladder the irreducible particle- 
particle scattering vertex Tj peaks at (tTjTt) momentum transfer as for the 2D Hubbard model. 
These results suggest that the peaking of F/ near (vr, tt) momentum transfer along with the enhanced 
single-particle spectral weight enhances the pairing correlations, and this is the reason for the strong 
dependence of D on t±/t seen in Fig 7.4. 



7.2.2 Irreducible particle-particle interaction 

In the previous sections, it has been seen that the 2-leg Hubbard ladder has power-law d^2_y2- 
wave pairing correlations in its ground state, and for a range of the model parameters the pairing 
correlations can get strong. It is useful to gain insight into the mechanism which leads to the 
d^2_y2-wave pairing correlations in this model. For this reason, here the Monte Carlo results on 
the irreducible particle-particle interaction F/ of the 2-leg Hubbard ladder will be reviewed. These 
are results from [Dahm and Scalapino 1997]. Comparisons with the magnetic susceptibility x(Q)O) 
will show that the short-range AF spin-fluctuations are responsible for the momentum structure in 
F/ for q near {n,7r). Using the data on F/ and the single-particle Green's function G, the Bethe- 
Salpeter equation will be solved in the particle-particle channel and the leading singlet pairing 
channel will be shown to have c/^2_„2-like symmetry. 
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(a) 




Figure 7.7: Distribution of the single-particle spectral weight A{k, to) in the k and lo plane. The 
intensity of the shading indicates the amount of the spectral weight. These results are for U = 4t, 
{n) = 0.94, T = 0.125t and (a) t±/t = 1.4, (b) 1.6 and (c) 1.8. 



Using Monte Carlo simulations, the singlet irreducible vertex Tis{fl,iiOm) has been calculated 
on a 2 X 16 lattice for U /t = A and (n) = 0.875 in the same way as discussed in Section 6.2. In 
addition, here t±/t = 1.5 was chosen so that the system has strong pairing correlations in the 
ground state. In Fig. 7.8(a), Vjs{f\,iijJm = 0) is plotted as a function of q where q = (q',7r) at 
various temperatures. Here, as in Section 6, q = p — p' is the momentum transfer and p' is kept 
fixed at ('7r,0) while p is scanned. At high temperatures, r75(q) is flat in momentum space with a 
magnitude varying between 8t and lOt, and as T is lowered to 0.25t, r/s(q, 0) develops significant 
amount of weight at q = (tTjTt) momentum transfer becoming of order 20t. This behaviour is 
similar to what has been seen in Section 6 for the 2D Hubbard model. 

Figure 7.8(b) shows the momentum dependence of the magnetic susceptibility xi^i^) foi" the 
2-leg Hubbard ladder for the same model parameters. Here, xi^^ 0) is also plotted as a function of 
q where q = (g, vr). Comparing with Fig. 7.8(a), one observes that the evolution of Tjs with tem- 
perature is closely related to that of x(Q)O)) which implies that the short-range antiferromagnetic 
correlations are responsible for the momentum structure in r75(q, 0) for q near (tt, tt). 
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Figure 7.8: Monte Carlo results on the momentum dependence of (a) the irreducible particle- 
particle interaction in the singlet channel Tis{ci,iuj„i = 0) and (b) the magnetic susceptibility 
x{<i,iu}m = 0). These results were obtained on a 2 x 16 lattice for U = At, T = 0.25t, {n) = 0.87 
and t± = 1.5t. 
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Figure 7.9: The leading singlet eigenfunction (f){p, iu)n) 
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Figure 7.10: (a) Temperature and (b) tj^/t dependence of the leading singlet eigenvalue Ai. These 
results were obtained on a 2 x 16 lattice for U = At and (n) = 0.87. In (a) t±/t = 1.5 and in (b) 
T = 0.25i were used. 

7.2.3 Bethe-Salpeter equation 

In order to see what type of pairing state is favored by the irreducible particle-particle vertex seen 
in the previous section, here results from the solution of the Bethe-Salpeter equation 

T 

P' 

from Ref. [Dahm and Scalapino 1997] will be shown. Figure 7.9 shows the momentum dependence 
of the leading eigenfunction 0(p, at tj„ = vrT for T = 0.25t. Figures 7.10(a) and (b) show the 
dependence of the leading singlet eigenvalue Ai on T/t and t^/t, respectively. Note that (j){p,i'irT) 
peaks near (7r,0) and (0,7r), and it changes sign between these two points. Hence, in this sense it is 
d^2 _y2-\ike, but it does not have a node since ^(p,z7rT) docs not vanish near any of the four Fermi 
surface points. In the previous section, wc have seen that for t±/t ~ 1.5. the bonding band has 
spectral weight near the Fermi level for p ~ (vr,0), and the antibonding band has spectral weight 
near the Fermi level for p ~ (0,7r). Hence, these Fermi points can be connected by scatterings 
involving q = (vr, vr) momentum transfer. Since Tjg is large and repulsive for q ~ (vr,7r), the 
leading singlet eigenfunction (p of the Bethe-Salpeter equation has opposite signs for p near (7r,0) 
and (0, tt). 

In Fig. 7.10(a), it is seen that the eigenvalue of the dr^2_y2-\ike eigenfunction increases by about 
an order of magnitude as T is lowered from It to 0.25t. However, even at T = 0.25t, which is the 
lowest temperature where the Monte Carlo calculation of the Bethe-Salpeter eigenvalues can be 
carried out, the leading singlet eigenvalue is only about 0.3. 

Figure 7.10(b) shows that at T = 0.25t the leading singlet eigenvalue has in fact a weak de- 
pendence on t±/t, and only a broad peak at t±/t = 1.6 is seen. This is unexpected considering 
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the strong dependence of D on t±/t seen in the previous section. It means that at T = 0.25t the 
pairing interaction has not reached its full strength and, in addition, the thermal smearing effects 
significantly weaken the pairing correlations which are observed in the ground state by the DMRG 
technique. 

7.3 Comparison of the 2-leg and the 2D Hubbard models 

It is interesting to compare these results on the 2-leg ladder with those on the 2D case seen in 
Section 6. In Fig. 7.10(b), one sees that at T = 0.25t, U/t = 4 and t±/t = 1.0 the leading singlet 
eigenvalue is 0.22, which is slightly larger than = 0.18 at the same U/t, (n) and T/t values for the 
2D case. For t±/t = 1.0, the DMRG calculations find that D{{.) varies approximately as in the 
ground state and, hence, the pairing correlations are only as strong as in the U = case. If one just 
uses this comparison of the eigenvalues for the 2-leg and the 2D cases, then one would expect weak 
pairing correlations in the ground state of the 2D Hubbard model. However, in 2D, for U = 8t 
is found to be higher than for U = At at T = 0.5t. Furthermore, from the DMRG calculations one 
knows that when t±/t is tuned to the right value, in the 2-leg ladder strong pairing correlations 
can develop. A similar dependence on the model parameters, such as the second-nearest-neighbour 
hopping t' , could exist in the 2D Hubbard model. 

It has been seen above that the QMC results on the irreducible particle-particle vertex of the 
2-leg ladder are similar to those on the 2D Hubbard model. In the ground state of the 2-leg ladder, 
the DMRG calculations find enhanced d.j.2 _y2-wave pairing correlations in a certain range of the 
parameters. Both in the 2D and the 2-leg models, we have seen that, when doped with holes there 
are short-range AF correlations and that they strongly influence the low-energy single-particle 
properties. In both cases, short-range AF correlations cause r7(q, 0) to peak at (7r,7r) momentum 
transfer. It needs to be noted that the pairing correlations observed in these two models do not 
require a particularly sharp peak in r7(q, 0) at q = (tt, tt), but rather simply weight at large 
momentum transfers. 

Currently, it is not possible to determine whether the doped 2D Hubbard model has long-range 
dx2_y2-wavc superconducting order in its ground state, or, if it did, whether it would be sufficient to 
explain superconducting transition temperatures as high as those found in the cuprates. In spite of 
this, the results discussed above show that effects which increase the single-particle spectral weight 
near the (tt, 0) and (0, tt) points of the Brillouin zone as well as effects which increase the strength 
of the particle-particle interaction at (tt, tt) momentum transfer will act to enhance dj.2_y2-wave 
pairing. 

It is possible that the strength of the dx2_y2 pairing correlations in the ground state of the 2D 
Hubbard model depends sensitively on the model parameters such as the second-near-neighbour 
hopping t' in a way similar to what is seen for the 2-leg ladder in Section 7. In the 2-leg case, the 
pairing correlations are as weak as those of the U = system when t± = t. However, when t±/t is 
tuned so that the fiat bands are located near the Fermi level, the system exhibits enhanced pairing 
correlations. It is possible that the 2D Hubbard model with only near-neighbour hopping similarly 
exhibits weak pairing correlations in the ground state, but there can be enhanced pairing when an 
additional parameter such as t' is tuned. In the 2-leg ladder case, the value of t±/t for which the 
fiat bands are near the Fermi level is renormalized by the Coulomb repulsion. Similarly, in the 2D 
case, the optimum value of t' could be renormalized, and it is difficult to estimate it in advance. In 
a bilayer Hubbard model, the bilayer coupling could also play a role. These are issues which need 
to be resolved by exact techniques in the future. 

Beyond these, one would expect that any additional contribution to the irreducible electron- 
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electron interaction which is repulsive for q ~ (vr, vr) momentum transfers or attractive for q 
would act to enhance the (1^2 _y2 pairing, when added to the 2D or the 2-leg Hubbard system. 
For instance, when a phonon mediated interaction which is most attractive for q ~ momentum 
transfers is added to an AF spin-fluctuation exchange interaction, it is found within the t-matrix 
approximation that the d^2_y2 eigenvalue gets enhanced [Bulut and Scalapino 1996]. However, it 
is necessary to study such effects using exact techniques. 

It could also be possible to design flat band dispersion near the Fermi level by controlling the 
lattice geometry and parameters [Imada and Kohno 2000]. Imada and Kohno have carried out 
exact diagonalization calculations for a ID 16-site t-J model with additional three-site terms and 
longer range hoppings. By tuning the longer-range hopping parameters, they have created fiat band 
dispersion near the Fermi level and, in this case, they find an enhanced spin gap and an enhanced 
tendency for pairing. They have also proposed various multiband models which could exhibit flat 
bands near the Fermi level and enhanced pairing. 

Even though for the 2D case the low-doping and the low-temperature regime where a spin gap 
could exist is beyond the reach of the exact techniques, in the 2-leg case the spin gap can attain 
large values. In the doped 2-leg ladder, is maximum when the bottom of the antibonding band 
is near the Fermi level [Noack et al. 1996]. For instance, for U = 8t, (n) = 0.875 and t± = IM, 
the spin gap has the value of 0.06t, which corresponds to ~ 300K for a t of order 0.45 eV. 

It is also useful to compare the density correlations seen in the 2-leg and the 2D Hubbard 
models. With the DMRG method [Noack et al. 1996], the following density-density correlation 
function has been calculated for the 2 x 32 Hubbard ladder, 

S{i,j,\) = {riixfiji) - {nix){nji). (7.8) 

Here, riix is the electron occupation number at the z'th site of chain A. By Fourier transforming, 
iS'(q) has been obtained for U = 8t, (n) = 0.875 and various values of t±/t. No obvious feature is 
found in S'(q) at the "2kF" wave vector of the 2-leg Hubbard ladder, which is q* = {q*,7r) with 
q* = 7r(n). On the other hand, a feature is observed at the "4kf " wave vector, which corresponds to 
(7r/4, 0) for (n) = 0.875. Especially for t±/t = 1.5, this feature becomes more obvious. In order to 
isolate the "4ki?" component of the density correlations, a correlation function A^(q) involving four 
density operators has been calculated. This correlation function exhibits a clear peak at (7r/4, 0) 
for both t^/t = 1.0 and 1.5, and in real space it decays as power law, while S{i,j,X) decays 
exponentially. 

The results on the 2-leg Hubbard ladder show that, when doped, this model exhibits simultane- 
ously short-range AF correlations and power-law decaying d^2_y2-\ike superconducting and 'Mk^" 
density correlations. For t±/t ~ 1.5, the superconducting correlations decay more slowly than 
the "4ki7'" density correlations. The 2D Hubbard model also exhibits short-range AF and dx2_y2 
superconducting correlations. In addition, the QMC data on the 2D case seen in Section 4 imply 
that the features found in the density susceptibility n(q, zw^ = 0) might be related to the 'Mk^" 
wave vectors rather than "2ki;'" for large U/t. Hence, the 2-leg and the 2D Hubbard models, when 
doped, appear to have similar magnetic, superconducting and density properties. Furthermore, 
in both cases, the single-particle dispersion near the (tt, 0) and (0, vr) points get flattened by the 
many-body effects. These flat bands also seem to play a key role in determining the strength of the 
pairing correlations in both models. At this point, it is necessary to note that in order to compare 
the 2-leg Hubbard model with the 2D case, one should not use isotropic hopping, t±/t = 1.0, 
since in this case the Fermi surface points are not connected by scatterings involving Q ~ (7r,7r) 
momentum transfers. Instead, anisotropic hopping, where the bonding and the antibonding Fermi 
surface points can be connected by (tt, tt) scatterings, need to be used. 
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Finally, an important question is whether the 2D Hubbard model exhibits spin-charge separation 
as in the ID case [Kivclson et al. 1987, Anderson et al. 1987]. While the 2D case cannot be resolved 
currently, in the 2-leg Hubbard ladder no indications of spin-charge separation as in the ID Hubbard 
model have been found [Noack et al. 1996]. 

8 Discussion 

The numerical results discussed above point out that the short-range AF spin fluctuations are re- 
sponsible for the dx2_y2-wme superconducting correlations. The relation of the AF spin fluctuations 
to the d^2_y2-w&Ye superconductivity has been studied using various diagrammatic approaches. 
Perhaps, the most commonly used of these approaches is the fluctuation exchange (FLEX) ap- 
proximation [Bickers et al. 1989], which self-consistently treats the fluctuations in the magnetic, 
density and the pairing channels. The FLEX technique has been used for obtaining possible phase 
diagrams and the estimates of the superconducting transition temperatures. In Section 8.1 below, a 
comparison of the QMC data with the results of the FLEX approach to the 2D Hubbard model will 
be carried out for the single-particle and the pairing properties. The purpose of this comparison is 
to have an idea for the range of applicability of the FLEX approximation, and see how it should 
be extended further. Another reason for carrying out such a comparison is because the Eliashberg 
type of calculations of the Tc's using the spin fluctuations for the cuprates [Monthoux and Pines 
1992] or the similar self-consistent spin-fluctuation exchange calculations [Moriya et al. 1990] are 
basically at the same level with the FLEX approximation. Hence, it is of interest to compare with 
the exact QMC calculations. 

It will be seen that for ?7 = 4f, the FLEX provides results in excellent agreement with the QMC 
data on the density of states N{u) and the Bethe-Salpeter eigenvalues in the pairing channel at the 
temperatures where the QMC calculations are carried out, as found earlier [Bickers et al. 1989]. 
For stronger coupling U = 8t, there are differences between the FLEX and the QMC results. The 
correlated metallic band which develops by doping the AF Mott-Hubbard insulator is not obtained 
within FLEX and this affects the strength of the dj.2_2^2-wave pairing correlations. 

Following this, in Section 8.2, other types of Monte Carlo approaches to the Hubbard model, in 
particular the variational and the projector Monte Carlo algorithms, will be briefly discussed. In 
Section 8.3, the results of the recent dynamical-cluster approximation and the one-loop RG calcu- 
lations for dy.2_y2 pairing in the 2D Hubbard model will be discussed. There is also much interest 
in understanding the unusual normal state properties of the cuprates in the pseudogap regime. 
The origin and the implications of the pseudogap are important issues in this field. In Section 8.4, 
various calculations on the low-doping regime of the 2D Hubbard model will be discussed briefly 
and their results will be compared with the pseudogap seen in the cuprates. 

The t-J model which is closely related to the Hubbard model has also drawn broad attention. 
It is useful to compare the QMC and the DMRG results on the Hubbard model with the numerical 
studies of the t-J model. Below, in Section 8.5, such a comparison will be given. Here, the 
attention will be concentrated on the density and the pairing properties since that is where there 
are unresolved issues. A question of interest is whether in the ground state of the doped 2D Hubbard 
model there are special density correlations as seen in the t-J model, for instance, phase separation 
or stripes, and, if so, whether they would suppress the d^2_y2-waye pairing correlations observed 
at higher temperatures with the QMC simulations. 

In Section 8.6, what these numerical studies of the Hubbard model imply for the mechanism of 
the dj.2_2^2-wave superconductivity seen in the high-Tc cuprates will be discussed. The issue which 
will be addressed here is whether the 2D Hubbard model, or some variation of it, is sufflcient for 
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explaining why the values of the Tc's found in the cuprates are so high. 

8.1 Comparisons with the fluctuation-exchange approach to the 2D Hubbard 
model 

The FLEX approximation was used first by Bickers et al. to study (ia.2_j^2-wave superconductivity 
in the 2D Hubbard model [Bickers et al. 1989, Bickers and White 1991]. This approach self- 
consistently incorporates the fluctuations in the magnetic, density and the pairing channels. It is 
an approximation around the band limit, and it is conserving in the sense that the microscopic 
conservation laws for the particle number, energy, and momentum are obeyed. Within FLEX, it 
is found that the (1^2 _y2 pairing correlations are mediated by the AF spin fluctuations. The phase 
diagram of the U = At Hubbard model within the FLEX approximation is shown in Fig. 8.1. The 
FLEX calculations find a finite mean-field Neel temperature at half-filling and at small dopings up 
to 6%. Neighboring the SDW phase is a dj.2_y2-wskYe superconducting phase which is stabilized for 
dopings between 6% and 20%. Note that in the 2D Hubbard model, a finite Neel temperature is 
not possible and only a Kosterlitz-Thouless type of superconducting transition can occur. Hence, 
the finite transition temperatures seen in Fig. 8.1 represent transitions to regimes where the d,^2_y2 
correlations have a power law decay at finite T. In order to induce true long-range order at finite 
temperature, three dimensional couplings would be required. 




Figure 8.1: Sketch of the phase diagram of the 2D Hubbard model within the FLEX approximation 
for U = At. From [Bickers et al. 1989, Bickers and White 1991]. 

In addition, using the FLEX approximation the effects of the nearest-neighbor hopping t' has 
been investigated for the 2D Hubbard model [Monthoux and Scalapino 1994, Dahm and Tewordt 
1995]. This approach has been also used for studying a three-band Cu02 model which has nearest 
neighbor copper-oxygen hopping and an onsite Coulomb repulsion at the Cu sites [Luo and Bickers 
1993]. The solutions of the FLEX equations in the dx2-y2-wsve superconducting state were also 
obtained [Pao and Bickers 1994 and 1995, Monthoux and Scalapino 1994, Dahm and Tewordt 1995]. 

In the following, a comparison of the FLEX results with the QMC data will be carried out using 



60 



results from Ref. [Dahm and Bulut 1996]. First the single-particle properties and then tlic pairing 
correlations will be discussed. Figure 8.2 compares the QMC data on the density of states N{uj) 
with the FLEX results for various values of U/t at (n) = 0.87. While for U = 4i, the results from 
the two different approaches are similar at these temperatures, there are qualitative differences for 
U = 8t and 12t. The Fermi level within FLEX moves slower with the doping at large U /t. The 
correlated metallic band at the Fermi level as well as the lower and the upper Hubbard bands and 
the Mott-Hubbard pseudogap are not obtained within the FLEX approximation. For [/ = 4t and 
T = O.SSt, the FLEX and the QMC results have similar qualitative features. However, at lower 
temperatures A^(uj) could develop a pseudogap for U = At also. Figure 8.3 shows the temperature 
evolution of N{lu) for U = 8t and (n) = 0.87 within the FLEX approximation. Comparing this 
with Fig. 5.3, one sees that the differences with the QMC data continue to exist as T is lowered. 
In Fig. 8.4, the QMC and the FLEX results on N{uj) are compared at half-filling for U = 8t and 
T = 0.5t. Here, it is seen that the development of the Mott-Hubbard gap at half-filling is not 
obtained within FLEX. 
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Figure 8.2: Comparison of the FLEX and the QMC results on N{uj) versus u at {n) = 0.875 for 

(a) U = At, (b) U = 8t, and (c) U = 12t. In (a) and (b) T = 0.33t was used and in (c) T = OM. 
Here, the vertical long-dashed and the short-dashed lines denote the chemical potential for the 
FLEX and the QMC calculations, respectively. 

Next, the FLEX results on the d^2_y2-wa,vc eigenvalue A^^ of the Bcthe-Salpeter equation in the 
particle-particle channel will be shown. In Fig. 8.5, is plotted as a function of T/t for U/t = 4 
and 8 in the temperature regime where the QMC data exist. Here, it is seen that as U/t increases 
from 4 to 8, A^ changes by a small amount. For U = 4t, the FLEX results on A^ versus T/t are in 
good agreement with the QMC data. However, for U = 8t the FLEX approximation underestimates 
\d by about a factor of two. This is the major difference between the QMC and the FLEX results 
on Arf. Within FLEX the effective pairing interaction is also attractive in the odd-frequency s and 
p-wave channels, in addition to the even-frequency d-wave channel. 

Comparing Fig. 5.3 with Fig. 8.3, one sees that the density of states at the Fermi level is also 
underestimated by about a factor of two within FLEX. From the simple spin-fluctuation exchange 
form of Eq. (2.1), one expects that the d^2_y2-wave pairing interaction increases with U/t. As 
discussed in Section 6.2, the QMC data also shows that the reducible vertex Fg gets enhanced as 
U/t is increased from 4 to 8, even though the irreducible vertex Tjg could not be obtained for 
this value of U/t. Hence, it must be that the FLEX approach underestimates the amount of the 
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Figure 8.3: Temperature evolution of N{iJ) versus uj within the FLEX approximation for U = 8t 
and (n) = 0.875. Here, the vertical long-dashed line denotes the chemical potential. 
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Figure 8.4: Comparison of the FLEX and the QMC results on N{u) versus u at half-filling for 
?7 = 8t and T = OM. 



single-particle spectral weight at the Fermi level, and this causes the d^2_y2-wave eigenvalue to be 
smaller with respect to the QMC value. 

Here, it is seen that the correlated metallic band which forms upon doping the AF Mott-Hubbard 
insulator is not obtained within FLEX, which appears to cause why the FLEX underestimates for 
U = 8t. Whether similar effects could take place when a CDW insulator is doped is an interesting 
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Figure 8.5: Temperature dependence of the d^2_y2-wave eigenvalue within the FLEX approxima- 
tion at (n) = 0.875 for U = At and 8t. 

problem. 
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Figure 8.6: Superconducting transition temperature Tc versus ?7/i within the FLEX approximation 
for (n) = 0.875. From [Pao and Bickers 1994]. 

Above, it has been seen that within FLEX has a weak dependence on U/t at temperatures 
between 0.25t and l.Ot. Related to this is Fig. 8.6 which shows the transition temperature Tc 
obtained within FLEX as a function of U/t. Here, it is seen that Tc has a broad peak at ~ 6t. 
In this section, it has been seen that the FLEX underestimates the strength of the d^2_y2-ws:ve 
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pairing when U is of order 8t. An important question raised by these comparisons is whether for 
U ^ 8t the Tc should be higher than 0.025t. There is a chance that the maximum possible Tc in the 
2D Hubbard model is higher than the FLEX estimate. However, note that these values for Tc are 
just estimates for a possible Kosterlitz-Thouless transition. In fact, the ground state of the doped 
2D Hubbard model is not known, and it is currently beyond the reach of the exact many-body 
calculations. 

8.2 Other Monte Carlo results on the 2D Hubbard model 

The QMC data presented in the previous sections were obtained with the determinantal Monte 
Carlo technique [White et al. 1989b]. Variational and projector Monte Carlo algorithms were also 
used for studying the 2D Hubbard model. However, these approaches have arrived at different con- 
clusions about the possibility of the dj.2_j^2-wave pairing. In this section, some of these calculations 
will be reviewed briefly. 

In order to see whether superconducting long-range order develops in the ground state of the 
doped 2D Hubbard model, zero temperature projector Monte Carlo calculations have been carried 
out [Imada 1991, Furukawa and Imada 1992]. In these calculations, an optimised initial state was 
used as a guiding function, mitigating the sign problem but making the calculation variational rather 
than exact. With this technique, the equal-time pair-field correlation functions were calculated for 
the s, extended s and the d^2_y2-wave singlet pairing channels. However, no size dependence was 
found in the data, which implies that there is no long-range superconducting order with these 
symmetries in the ground state. 

The constrained-path Monte Carlo (CPMC) algorithm developed in Ref. [Zhang et al. 1995] 
is a variational method which starts from a trial wave function and uses the same |^'t) as a 
constraining wave function as the simulation is carried out in the space of the Slater determinants. 
In calculations on up to 16 x 16 lattices, free-electron and unrestricted Hartree-Fock wave functions 
were used as j^'r); however no d^2_y2-wave superconducting long-range order was found [Zhang 
et al. 1997]. The CPMC method has been extended by using as the constraining function a 
d^2_y2-wave BCS wave function [Guerrero et al. 1999]. Superconducting long-range order was not 
found in this case either. An important issue is how well the constraining function describes the 
correlations in the Hubbard model. In fact, a comparison of the CPMC and the DMRG results, 
which will be discussed later in Section 8.5.1, has been carried out for the magnetic and the density 
correlations in the 12 x 3 Hubbard model [Bonca et al. 2000]. This comparison has found that the 
CPMC results depend sensitively on the constraining wave function especially for large U/t ^ 8. 
In this respect, it would be useful to test how well the CPMC method describes the dx2_y2-\ike 
superconducting correlations found in the 2-leg Hubbard ladder. For instance, it is known from the 
DMRG calculations that in the 2-leg Hubbard ladder for U/t = 8, t±/t ~ 1.5 and near half-filling 
the pair-field correlation function D{£) decays as , where < 1.0. It would be useful to see 
whether the CPMC method describes the power-law pairing correlations which exist in the 2-leg 
Hubbard ladder. 

Another variational Monte Carlo approach uses as the trial wave function [Nakanishi et al. 
1997, Yamaji et al. 1998] 

I^t) = PnPgI'^bcs), (8.1) 

where 

\^BCS) = + Vkcl^cl^^^)\0) (8.2) 
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is the usual dj.2_j^2-wave BCS wave function, 

PG = Ui[l-{l-g)n,^ni^], (8.3) 

is the Gutzwiher projection operator, and P/v is an operator which projects out states with electron 
number N. In this method, the magnitude of the d^2_y2-wave gap function entering the usual BCS 
coefficients Uk and Uk, the parameter g, and the chemical potential /x are used as variational 
parameters, and their values are determined by minimising the ground state energy with a Monte 
Carlo simulation. With this method it is found that, for fillings between 0.84 and 0.68, the dj.2_y2- 
wave state has the lowest energy, and for (n) > 0.84 an SDW state is favored. In addition, when 
a next-near-neighbor hopping term t' is turned on, the superconducting condensation energy gets 
enhanced. For instance, for U = 8t this enhancement is maximum when t' ~ — O.lt. It should 
be noted that this method was also applied to the 2- leg Hubbard ladder [Koike et al. 2000]. In 
this case, the gap function was assumed to take momentum-independent values Ai and A2 on 
the bonding and the antibonding bands. Treating Ai and A2 as variational parameters, it was 
found that a superconducting state with dx2_y2-\\\.e symmetry is stabilized when the bottom of 
the antibonding band is near the Fermi level, which is in agreement with the exact-diagonalization 
[Yamaji and Shimoi 1994] and the DMRG [Noack et al. 1995 and 1997] calculations. 

Next, a projector Monte Carlo approach [Husslein et al. 1996] will be discussed where the 
ground state wave function of the 2D Hubbard model is estimated from 

=e^^|*o). (8.4) 

Here, |^'o) is the ground state of the nonintcracting electrons and the parameter 9 is taken to be 
about 8 on a 12 X 12 lattice. This approach has been used to calculate the (i^2_j,2-wave pair-field 
correlation function in the weak-coupling regime, 0.5t <U < 3t, near half-filling. It is found that 
the system has long-range d^2_y2-wme pairing correlations for negative values of t' , for instance, 
at t' = -0.3t for U = 2t. 

In this section, it was seen that various approaches arrive at different conclusions about d^2_y2- 
wave pairing in the 2D Hubbard model. This emphasises the importance of the choice of the 
trial wave functions and points out at the need for carrying out calculations without introducing 
approximations. 

8.3 Dynamical cluster approximation and RG calculations for dx2^y2 pairing in 
the 2D Hubbard model 

Recently, the dynamical cluster approximation [Maier et al. 2000] and the one-loop renormalization- 
group method employing a 2D Fermi surface [Zanchi and Schulz 2000, Halboth and Mctzner 2000, 
Honerkamp et al. 2001] were used for studying dj.2_y2 pairing in the 2D Hubbard model. In this 
section, the findings of these studies will be discussed briefly. 

When the infinite dimensions limit of the Hubbard model is taken with proper scaling, the 
many-body problem becomes local [Metzner and Vollhardt 1989, Miillcr-Hartman 1989], and it can 
be mapped to an Anderson impurity problem, which can then be solved with various many-body 
techniques [Pruschke et al. 1995, Georges et al. 1996]. This is called the dynamical mean-field 
approximation (DMFA). The DMFA is interesting because it is a strong coupling technique. For 
instance, for large U/t and off of half-filling, the DMFA yields a narrow peak in N{uj) near the 
Fermi level [Jarrell 1992], which is also seen in the 2D QMC data. However, the DMFA does not 
incorporate the non-local correlations, and hence it is not possible to study d^2_y2 pairing with 
it. The DC A incorporates the non-local corrections to DMFA by mapping the lattice problem 
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onto an embedded cluster of size Nc, rather than onto an impurity problem. In DCA calculations, 
the dynamical correlation length is restricted to L = y/TT^, and the DCA would become exact 
for Nc ^ CO, while it reduces to DMFA for Nc = 1. With this approach, the mean-field d^2_y2 
superconducting Tc's were calculated for A^c = 4 by using the non-crossing approximation to solve 
the cluster problem [Maier et al. 2000]. For U = 12t, it is found that Tc has the maximum value of 
~ 0.05t for about 20% doping. It is also found that Tc increases for positive values of the second- 
nearest-neighbour hopping t' , and decreases for negative values of t' . This result agrees with the 
DMRG calculations on the t-J model which find d^2_y2 pairing for t' > [White and Scalapino 
1998a]. With DCA, the low doping regime of the 2D Hubbard model was also studied, and these 
results will be discussed in Section 8.4. 

After the discovery of the high Tc cuprates, the one-loop RG approach was extended from ID to 
2D in order to study this problem [Dzyaloshinskii 1987, Schulz 1987, Lederer et al. 1987]. The RG 
calculations are interesting, because with this technique the particle-particle and the particle-hole 
channels are treated on equal footing. These studies focused on the scattering processes between 
the Fermi surface regions near the van Hove singularities. For the 2D Hubbard model, the AF 
SDW state at half-filling and the (1^2 _y2 superconducting state induced by AF fluctuations away 
from half-filling were found [Schulz 1987, Lederer et al. 1987]. The scatterings involving the full 
2D Fermi surface were taken into account with the work of Zanchi and Schulz, who studied the RG 
flows using a 32-patch discretization of the 2D Fermi surface for = [Zanchi and Schulz 2000]. 
These calculations found two different regimes, one dominated by the AF correlations and the other 
by the (1^2 _y2 pairing. The 2D RG calculations were also extended to the t' 7^ case [Halboth and 
Metzncr 2000]. In the calculations for t' 7^ by Honerkamp et al., an intermediate regime is 
found between the two regimes dominated by the AF correlations and d^2_y2 pairing [Honerkamp 
et al. 2001]. In this intermediate regime, there are competing AF and d^2_y2 superconducting 
correlations. This regime exists only for sizeable t' < and it exhibits features which are similar to 
those seen in the pseudogap regime of the cuprates, which will be discussed in the following section. 
Honerkamp et al. extracted a temperature versus doping phase diagram from the 2D RG flows for 
the 2D Hubbard model with t' < 0, which is similar to that of the cuprates. 

In spite of these interesting results it has to be kept in mind that the dynamical correlation 
length in DCA is cut-off by the size of the cluster and the one-loop RG is a weak-coupling approach. 
In addition, in RG calculations the single-particle self-energy corrections are not included at the 
one-loop level. In Section 8.1, it was seen that the single-particle self-energy corrections could play 
an important role in determining the strength of pairing. It should also be noted that, while the 
RG finds that the regime with t' < is favored, in the DCA calculations the mean-field Tc's are 
higher for t' > 0. These are some of the reasons for why exact results are necessary in order to 
reach a final conclusion. 

8.4 Low-doping regime of the 2D Hubbard model 

In the normal state of the underdoped cuprates, a pseudogap is seen in the excitation spectrum 
below a temperature T* which depends on the doping. The nature of the pseudogap is an important 
problem in this field and there exist a wide range of ideas about its origin and its implications 
[M^S-HTSC VI proceedings 2000]. In this section, various calculations which have been carried 
out for exploring whether there is a pseudogap regime in the 2D Hubbard model will be discussed. 
In Ref. [Jaklic and Prelovsek 2000], a review of the numerical calculations on the t-J model at 
finite temperatures is given and the results of these calculations are compared with the anomalous 
normal-state properties of the cuprates including the pseudogap. 
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The normal state pseudogap, while exhibiting dependence on the material properties, is seen in 
the uniform magnetic susceptibility, the low-frequency optical conductivity, the ARPES spectrum 
and various other measurements of the electronic properties. An interesting feature of the normal 
state pseudogap observed in the ARPES spectrum is that it is anisotropic on the Fermi surface: 
the pseudogap has maximum amplitude near the (vr, 0) and (0, vr) points of the Brillouin zone and 
it is minimum for wave vectors near (7r/2,7r/2) [Ding et al. 1996, Ronning et al. 1998]. 

Kampf and Schrieffer first discussed the possibility of a pseudogap within the spin-bag approach 
for the 2D Hubbard model [Kampf and Schrieffer 1990a and 1990b]. They showed how a pseudogap 
could develop because of the AF fluctuations already at the one-loop level for the single-particle 
self-energy as the AF instability is approached. In the FLEX calculations, it was also found that 
a pseudogap opens in the density of states when the strength of the short-range AF fluctuations 
increases [Dahm and Tewordt 1995]. 

In spite of the sign problem, there are QMC data on A(p,u;) at temperatures as low as 0.25t in 
the underdoped regime [Preuss et al. 1997]. These data show indications that a pseudogap opens 
in A{p,u;) as T decreases. In particular, for (n) ~ 0.95, it is seen that spectral weight above the 
Fermi level for p between (vr, 0) and (tt, vr) decreases gradually as T decreases from OM to 0.25t. 
In these calculations, the pseudogap has been attributed to the AF spin correlations, which are 
becoming larger than a lattice spacing for T < 0.3t in the underdoped regime. However, in order 
to be able to make direct comparisons with the experimental data on the pseudogap, it would be 
necessary to reach temperatures below O.lt. 

In recent one-loop RG calculations for the Hubbard model with t' < 0, a saddle-point regime is 
found between the AF ordered and the d,j.2_y2 superconducting regimes [Honerkamp et al. 2001]. 
This regime exhibits features similar to those seen in the pseudogap regime of the cuprates. Here, 
the uniform magnetic and charge susceptibilities flow to zero because of the pairing and the 
umklapp-scattering processes, respectively. The possible existence of an umklapp-gapped spin- 
liquid phase was suggested in earlier RG calculations where a two-patch discretization of the Bril- 
louin zone was used [Furukawa et al. 1998]. The similarity of this regime to the spin-gapped phase 
in the 2-leg Hubbard ladder was also noted. Because of the umklapp scatterings, in this regime, the 
sections of the Fermi surface which are near the (vr, 0) and (0, vr) points are truncated while at wave 
vectors near (vr/2,vr/2) gapless single-particle excitations remain. Clearly, these results are useful 
for interpreting a number of experimental data on the underdoped cuprates and, in particular, the 
ARPES data. 

The DCA calculations have also found interesting results about this subject [Huscroft et al. 
2001]. Using this method, the low-doping and the low-temperature regime of the 2D Hubbard 
model was studied. Here, QMC simulations were used to solve the embedded lattice problem 
with Nc = 8. It was found that, as T is lowered, an anisotropic pseudogap develops in A(p,Lj) 
and, in addition, the uniform magnetic susceptibility x(r) gets suppressed. In particular, for 
(n) = 0.95, U = 6t and T = 0.06t, a pseudogap is found in A(p,u!) for p = (vr, 0) and (0,vr) and 
not for (vr/2, vr/2). For these values oiU/t and (n), a mean-field d^2_y2 superconducting transition 
temperature of Tc ~ 0.04t is obtained with the DCA. In addition, a sharp drop in x{T) is observed 
near the temperature where the pseudogap in A{p, u) opens. In these calculations, the downturn 
in x{T) has been attribTitcd to the development of the AF correlations. 

These results support the idea that perhaps the normal state pseudogap seen in the cuprates 
could be understood within a 2D Hubbard framework. However, there still exist a wide range of 
ideas about the origin of the pseudogap, and it is one of the important unresolved issues in this 
field. 
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8.5 Comparisons with the t-J model 

In the large U/t limit, the Hubbard model reduces to [Hirsch 1985] 

H=-t (4,<TCj,a + c]^„Ci^a) + J^(Si • Sj - ^niUj) 

i,u 5^5' 

where J = At'^/U, and the double occupancy of a site is not allowed. In Eq. ( |8.5| ), i + 6 and i + 6' 
sum over the nearest neighbors of site i. The last sum in this expression, which involves operators 
acting at three different sites, is dropped and what is remaining is called the t-J model. Clearly, 
the t-J and the Hubbard models have differences. The numerical studies of the t-J model have 
produced valuable information about the magnetic, density and the superconducting correlations 
in this system. Reviews of these studies can be found in Refs. [Dagotto 1994, Jaklic and Prelovsek 
2000]. In this section, the density and the pairing correlations in the t-J and the Hubbard models 
will be compared. 

In Section 8.5.1, the results on the phase separation and the density correlations in the 2D t-J 
model obtained with various numerical techniques will be discussed. These data will be compared 
with the CPMC and the DMRG calculations for the 12 x 3 Hubbard lattice. In Section 8.5.2, the 
nature of the pairing correlations seen in the 2D t-J model will be compared with those in the 2D 
Hubbard model. The results on the 2-leg t-J and Hubbard ladders will be compared in Section 
8.5.3. 

8.5.1 Comparisons with the density correlations in the 2D t-J model 

The QMC simulations [Moreo and Scalapino 1991] and the exact-diagonalization [Dagotto et al. 
1992b] calculations on the 2D Hubbard model did not find any indication of phase separation of 
the system into hole-rich and hole-poor regions. The 2-leg Hubbard ladder does not show any 
evidence for phase separation either [Noack et al. 1994]. In contrast to the Hubbard model, the 
2D t-J model phase separates at any electron filling for sufficiently large values of the interaction 
strength J/t. Various techniques including the variational and the exact-diagonalization calcula- 
tions [Emery et al. 1990], the high-temperature series expansions [Putikka et al. 1992], and the 
exact-diagonalization calculations [Dagotto 1994] were used for studying phase separation in the 
t-J model. Recently, the Green's function Monte Carlo (GFMC) [Hellberg and Manousakis 1997, 
1999 and 2000, Calandra et al. 1998] and the DMRG [White and Scalapino 1998a, Rommer et al. 
2000] were used for obtaining the phase-separation boundary in the t-J model. Currently, where 
the true phase-separation boundary lies is controversial, since the approaches noted above arrive 
at different conclusions. 

Emery et al. suggested that the 2D t-J model phase separates at all interaction strengths 
[Emery et al. 1990]. Subsequent calculations found that the phase separation occurs only for 
J > t, as indicated by the solid curve in Fig. 8.7 [Putikka et al. 1992]. However, recently, the 
GFMC calculations [Calandra et al. 1998] suggested that the phase-separation boundary might 
occur at lower J/t values near half-filling. This is indicated by the empty circles in Fig. 8.7. In 
these calculations, next to the phase-separation boundary lies the regime where the doped holes 
form dx2_y2-wa,ve pairs. The recent GFMC calculations by Hellberg and Manousakis, on the other 
hand, find that the critical value of J/t for phase separation extrapolates to zero at low dopings as 
shown by the dashed curve in Fig. 8.7 [Hellberg and Manousakis 1997, 1999 and 2000]. In these 
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Figure 8.7: Sketch of the phase-separation boundary for the 2D t-J model obtained from vari- 
ous calculations. The solid line represents the phase-separation boundary obtained by the high- 
temperature series expansion [Putikka et al. 1992]. The open circles represent the results of the 
GFMC calculations from Ref. [Calandra et al. 1998]. The dashed line is from Ref. [Hellberg and 
Manousakis 1997], which was also obtained by the GFMC. The arrow indicates the DMRG result 
for the critical value of J/t for phase separation on a six-leg t-J ladder in the limit of zero doping 
[Rommer et al. 2000]. 

calculations, the t-J model phase separates in the parameter regime appropriate for the cupratcs. 
For instance, at 15% hole doping the critical value of J/t is about 0.4. The DMRG calculations 
find that, in the physically relevant regime, the ground state of the t-J model is striped and not 
phase separated [White and Scalapino 1998a, 2000]. The DMRG calculations are carried out on 
large lattices with periodic boundary conditions at the long edges and open boundary conditions 
at the short edges. In these calculations the phase separation occurs for J > t. For instance, in 
the six-leg Hubbard ladder, the critical value of J/t is ~ 1.4 as the doping approaches zero, which 
is indicated by the arrow in Fig. 8.7. 

The stripes observed in the DMRG calculations are a domain wall of holes across which there is 
a TT-phase shift in the AF background. According to the DMRG calculations, the stripe formation 
represents an instability of the system where pairs of bound holes combine to form domain walls. In 
particular, through the stripe correlations the system reduces the frustration of the AF background 
and lowers the kinetic energy of the holes. However, the issue of the phase separation versus 
the stripe formation in the 2D t-J model remains controversial. The role of the open boundary 
conditions in producing the static stripes was questioned [Hellberg and Manousakis 2000] , and the 
need to have more than one pair of holes in order to see the stripe correlations was noted [White and 
Scalapino 2000]. The calculation of the dynamical spin and charge susceptibilities for sufficiently 
large t-J systems would resolve these issues. 

The possibility of stripe formation in strongly correlated systems were noted in mean-field cal- 
culations of the 2D Hubbard model soon after the discovery of the cuprates [Zaanen and Gunnarson 
1989, Poilblanc and Rice 1989, Schulz 1989, Machida 1989]. There has been a surge of interest 
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in this field since the observation of static stripe ordering in the neutron scattering experiments 
on Lai.6-xNdo.4Sra;Cu04 [Tranquada et al. 1995]. Beyond mean-field, the nature of the stripe 
correlations in the Hubbard model were studied with the DMRG and the CPMC techniques for 
the 3-leg Hubbard ladder [Bonca et al. 2000] . The findings of these calculations are important and 
they will be briefly described here. 

The CPMC is an approximate method which projects the ground state from a trial state as has 
been discussed in Section 8.2. Both the DMRG and the CPMC calculations find that for U > 6t 
there are static stripes in the ground state of a 12 x 3 Hubbard lattice doped with six holes when 
open boundary conditions are used. In Figures 8.8 and 8.9, the rung magnetization 

S^{i) = Y,{S'ihj)) (8-6) 

and the rung hole density 

pii) = i2(P('^3)) (8.7) 

are plotted as a function of the rung location i at U = 8t. Here, S^{i,j) and p{i,j) are the spin and 
hole-density operators at the j'th site of the i'th rung. In Fig. 8.8, it is seen that between rungs 
3 and 4 the spins are ferromagnetically coupled, causing a vr-phase shift, and the magnetisation 
density is small. The same occurs at rungs 9 and 10. In Fig. 8.9, it is seen that at these sites 
the holes form domain walls. In these figures the results of the unrestricted Hartree-Fock (UHF) 
calculations are also shown. The UHF approximation produces results in agreement with the CPMC 
and the DMRG calculations when the reduced value of [/ = 3t is used. This behavior is similar 
in spirit to using a reduced Coulomb repulsion within RPA for the spin susceptibility as discussed 
in Section 3.3. It is important to note that the structure of the stripes seen in these calculations 
for U > 6t is quite similar to those found in the three- leg t-J ladder [White and Scalapino 1998b]. 
However, for U < 6t these features disappear and no evidence is found for static stripes in the 
ground state. In this regime, only some evidence for low-lying states with stripes are found. At 
weak U/t, the density of virtual holes due to the double occupancy increases, and, it has been 
noted that this might weaken the stripe correlations. At this point, it is of interest to explore 
whether there is a relation between the stripe patterns which are observed within the presence of 
the open boundary conditions and any possible "AIlf" CDW correlations. It should also be noted 
that the CPMC method does not show evidence for stripe correlations even at large U/t when 
the free-electron trial wave function is used instead of the UHF wave function. Hence, caution is 
necessary in choosing a trial wave function which is suitable for the ground state. 

These results were found using open boundary conditions. The DMRG calculations are not 
carried out with periodic boundary conditions, but the CPMC calculations can be. When carried 
out with the periodic boundary conditions on square lattices, the CPMC calculations do not find 
stripes in the hole and spin densities and in the hole and spin correlation functions. This indicates 
that the static stripe pattern seen in the ground state results from the open boundary conditions 
breaking the translational symmetry. With the periodic boundary conditions, there might be 
low-lying states with stripe patterns but it might be difficult to detect them numerically. Hence, 
whether static stripes exist in the ground state of the 2D Hubbard model with periodic boundary 
conditions remains as an important open question. Here and in Section 4, it has been seen that 
there are similarities in the density correlations of the t-J model and the Hubbard model with large 
U/t. But, there are differences as well. 
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Figure 8.8: Rung-magnetisation versus the rung location for the 12 x 3 Hubbard lattice doped 
with six holes. Prom [Bonca et al. 2000]. 




Figure 8.9: Rung-hole density versus the rung location for the 12 x 3 Hubbard lattice doped with 
six holes. Prom [Bonca et al. 2000]. 
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8.5.2 Comparisons with the superconducting correlations in the 2D t-J model 

When two-holes are doped into the half-filled t-J model, they form a dj.2_j^2-wave bound pair for 
interaction strengths relevant to the high-T^ cuprates, J/t ^ 0.35 [Poilblanc 1993, Dagotto 1994]. 
However, much controversy exists over what happens when more than one pair of holes is doped. 
The GFMC calculations by Calandra et al. find that they form a dr^2_y2-wmQ superconducting 
ground state for J jt ^ 0.35 and dopings relevant to the high-Tc cuprates [Calandra et al. 1998]. 
However, the GFMC calculations by Hellberg and Manousakis find that this parameter regime lies 
right at the boundary for the phase separation of holes [Hcllbcrg and Manousakis 1997 and 1999]. 

The DMRG calculations find that in the same regime the ground state has static stripes, and the 
system does not have (i^2_j^2-wave superconducting long-range order [White and Scalapino 1998a, 
1999]. The fact that the static stripes compete with the superconductivity is an important result of 
these calculations. However, when the static stripe patterns are broken by a second-near-neighbor 
hopping t' > 0, it is seen that the system develops long-range d^2_y2-wave superconducting order, 
and in this regime the low-lying stripe correlations coexist with the long-range superconducting 
order. 

Both the t-J and the Hubbard models exhibit d,j.2_y2-wave pairing correlations due to the AF 
spin fiuctuations. However, there are differences between these models, one of them being the no- 
double-occupancy constraint in the t-J model. It is possible that this constraint is more favorable 
of the stripe correlations and, because of this, the one-band Hubbard or a three-band Cu02 model 
might have weaker tendency for stripe formation than the t-J model [Daul et al. 2000, Jeckelmann 
et al. 1998]. 

Whether the doped 2D Hubbard model develops static stripe patterns in the ground state is 
an important issue. In the previous section, we have seen from the QMC data that the doped 2D 
Hubbard model has d^2_y2-wave pairing correlations which develop as the temperature is lowered 
and they are not weak correlations. However, if indeed a static striped phase is more favored 
in the ground state then this could suppress the growth of the superconducting correlations at a 
temperature lower than where the QMC simulations arc carried out. In this respect, the comparison 
of the CPMC and the DMRG calculations for the 12 x 3 Hubbard lattice are useful [Bonca et al. 
2000]. However, further study is necessary before reaching conclusions about the ground state of 
the doped 2D Hubbard model. 

The results shown here emphasize the interplay of the magnetic, density and the pairing corre- 
lations in the Hubbard and t-J models. In the cuprates, the substitution of nonmagnetic impurities 
gives useful information about the interplay of these correlations in these materials. It is known 
that, when substituted in place of planar Cu, nonmagnetic impurities strongly suppress the super- 
conducting Tc [Xiao et al. 1990] and locally enhance the AF correlations [Mahajan et al. 1994 
and 2000]. To the extend that a nonmagnetic impurity can be considered as a pure potential scat- 
terer, these experimental data give information about the magnetic and pairing response of these 
materials when perturbed in the density channel. 

8.5.3 Comparisons with the 2-leg t-J ladder 

The Lanczos calculations on the 2-leg t-J ladder showed that this system, which has a spin gap in 
the insulating state, can exhibit superconducting correlations upon doping [Dagotto et al. 1992a]. 
The mean-field calculations noted the d^2_y2-like structure of the superconducting pairs in the t-J 
ladder [Gopalan et al. 1994]. The DMRG calculations on long t-J ladders established that in the 
ground state there are power-law decaying superconducting correlations [Hayward et al. 1995]. In 
Ref. [Dagotto and Rice 1996], the properties of n-leg spin-1/2 ladders are reviewed and comparisons 
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are made with various ladder compounds. 

There are differences in the nature of the d.j.2_y2-like superconducting correlations seen in the t-J 
and the Hubbard ladders. For instance, in the Hubbard ladder for t±/t ^ 1.0, the superconducting 
correlations are weak, and the pair-field correlation function D{1) decays as £~^, which is like the 
U = case [Noack et al. 1994]. Only when the distribution of the single-particle spectral weight 
is such that it can make use of the momentum structure in F/, does the system exhibit strong 
pairing correlations. On the other hand, in the isotropic t-J ladder, the system has strong pairing 
correlations. For instance, D{1) decays approximately as t'^ for (n) = 0.8 and J/t = 1 [Hayward 
et al. 1995]. In addition, Schulz has shown that D{t) decays as in the limit (n) — ^ 1 [Schulz 

1999]. Hence, for isotropic hopping the t-J ladder has stronger pairing correlations compared to 
the Hubbard ladder. Probably, this is related to the fact that in the t-J model the exchange term 
J, which is an effective attractive interaction between the nearest-neighbor antiparallel spins, is 
introduced by hand. On the other hand, in the Hubbard model, the effective attractive interaction 
which is responsible for the pairing is generated by the onsitc Coulomb repulsion through higher- 
order many-body processes only when the system has the suitable conditions. 

8.6 Implications for (ia.2_j^2-wave pairing in the cuprates 

The presence of (1^2 _y2 pairing correlations in the Hubbard model and the nature of the effective 
interaction mediating it are issues which are of interest for studying high temperature supercon- 
ductivity in the cuprates. This is so especially after it became clear that the superconducting gap 
function in the high Tc cuprates is of the d^2_j^2-wave type. Within this context, enormous amount 
of research has been carried out on this model. 

In this article, the numerical studies of dj.2_j^2-wave pairing in the Hubbard model have been 
reviewed. The Monte Carlo simulations have shown the presence of the d^2_y2-wa,ve pairing corre- 
lations in the Hubbard model, even though sufficiently low temperatures, where long-range order 
could establish, have not been reached. For U/t = 4 and at the temperatures where the simulations 
are carried out, the values of the d^2_y2-wsve eigenvalues of the particle-particle Bethe-Salpeter 
equation are in agreement with the FLEX calculations. When carried out at low temperatures, 
the FLEX calculations find T^s of order 0.025i. Since t is estimated to be about 0.45 eV for a 
single-band model of the cuprates [Hybertson et al. 1990], this value of the Tc corresponds to 
~ 130 K. This is a high value reflecting the electronic energy scales of the model. Furthermore, 
the comparison of the FLEX and the QMC data suggests that the pairing could be stronger for 
U = 8t. However, these are only mean-field estimates for a possible Kosterlitz-Thouless transition 
and, in fact, the ground state of the doped 2D Hubbard model is not known. It might be that 
an additional parameter needs to be tuned in order to create optimum conditions for d^2_y2-wave 
superconductivity in two dimensions. For instance, in the 2-leg ladder case, it was necessary to tune 
t±/t. In this respect, the second-nearest-neighbour hopping t' may play a role for the enhancement 
of the dx2_y2 pairing in the 2D system like t± does for the 2-leg ladder. It is also possible that a 
three-band Cu02 model with an onsite Coulomb repulsion at the Cu sites offers additional degrees 
of freedom for creating more favourable conditions for pairing. 

The QMC results reviewed here show that the short-range AF correlations are responsible for the 
dx2_y2 pairing correlations in the 2D Hubbard model. These results also suggest that effects which 
enhance the magnitude of r/s(q ~ (tt, tt), = 0) and the single-particle spectral weight near the 
(vr, 0) and the (0,7r) points of the Brillouin zone act to enhance the d^2_y2-waYC superconducting 
correlations. Clearly, other effects which enhance dj.2_j^2-wave pairing can exist. Exploring these 
possibilities is an active research field. 
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9 Summary and conclusions 



In this paper, the numerical studies of the 2D and the 2-leg Hubbard models have been reviewed. 
For the 2D Hubbard model, data from the QMC simulations have been shown. These data represent 

what has been obtained over the years in the parameter regime allowed by the sign problem. For 
the 2-leg Hubbard ladder, the QMC results at finite temperatures and the DMRG data on the 
ground state have been shown. 

Here, the emphasis has been placed on the d^2_y2-wave superconducting correlations observed 
in the Hubbard model. In order to develop an understanding of the origin of these correlations, 
results on the spin, charge and the single-particle excitations have been shown along with the data 
on the particle-particle and the particle-hole interactions. The observation of the short-range AF 
fluctuations by the NMR and the neutron scattering experiments, and the unusual single-particle 
spectrum seen in the ARPES data are properties which support using a Hubbard framework for 
studying the pairing correlations of the cuprates. 

In the 2D Hubbard model, upon hole doping the long-range AF order is destroyed and the 
system exhibits short-range AF correlations. The maximum-entropy analysis of the QMC data 
shows that the AF and the Coulomb correlations strongly affect the single-particle properties. In 
particular, it has been seen that, as the strength of U/t increases, significant amount of single- 
particle weight remains pinned near the Fermi level at the (vr, 0) and (0, vr) points of the Brillouin 
zone. These generate phase space for scatterings in the d3;2_j,2-wave BCS channel. The results on 
the particle-particle and the particle-hole irreducible interactions have been also presented. It has 
been seen that, for U = 4t, a properly-renormalized single spin-fluctuation exchange interaction can 
describe the momentum, the Matsubara-frequency and the temperature dependence of the effective- 
particle-particle interaction. This means that, for these values of U/t and T/t, the attraction in 
the d^2_y2 channel is mediated by the AF spin fluctuations. 

The DMRG results on the 2-leg ladder, which have been presented in Section 7 are important. 
These calculations represent the first example where a purely repulsive onsite interaction leads to 
superconducting correlations which are enhanced over the noninteracting [U = 0) case in the ground 
state of a bulk system. For certain values of and in the intermediate coupling regime, this 
model exhibits enhanced power-law d^2_j^2-like superconducting correlations. The QMC simulations 
for this system show that the effective particle-particle interaction F/ peaks at (vr, vr) momentum 
transfer and, in addition, the d^2_y2 pairing correlations are strongest when there is significant 
amount of single-particle spectral weight near the Fermi level at the (tt, 0) and (0, vr) points in the 
Brillouin zone. This way the system makes optimum use of the momentum structure in F/ for 
d^2_y2-'Sffme pairing. 

These data were also compared with various other approaches to the Hubbard model, such as 
the diagrammatic FLEX approximation, the variational and the projector Monte Carlo simulations, 
the dynamical cluster approximation and the one-loop RG calculations. In addition, the similarities 
and the differences with the t-J model were briefly discussed. The implications for the dj.2_y2-wme 
superconductivity seen in the cuprates were also noted. 

It is difficult to reach conclusions about the strongly correlated systems. Nevertheless, here, the 
following conclusions are given: 

(1) An onsite Coulomb repulsion can lead to superconducting correlations which are enhanced 
with respect to the noninteracting {U = 0) case in the ground state of a bulk system. This is 
proven for the ease of the 2-leg Hubbard ladder. 

(2) The 2D Hubbard model exhibits d^2_y2-wa.ve pairing correlations which grow as the temper- 
ature is lowered in the parameter regime where the QMC simulations are carried out. The fastest 
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growing pairing correlations occur in the singlet dj.2_j^2-wave channel. These correlations are not 
weak. 

(3) The QMC simulations also find that the effective pairing interaction in the 2D and the 2-leg 
Hubbard models in the intermediate coupling regime and at temperatures greater or of order J/2 
is consistent with the spin-fluctuation exchange approximation. 

(4) Two factors which create optimum conditions for d^2_y2 pairing in the Hubbard model are 
enhanced single-particle spectral weight at the Fermi level near the (tt, 0) and (0, vr) points in the 
Brillouin zone and large weight in the effective pairing interaction F/ at large momentum transfers. 
Clearly, there can be other ways of enhancing the dr^2_y2 pairing correlations. 

Acknow^ledgment s 

Most of the numerical data presented here are from various calculations carried out with T. 
Dahm, R.M. Noack, D.J. Scalapino and S.R. White. The QMC data shown here were obtained at 
the San Diego Supercomputer Center. 

10 Appendix 

The results shown in Sections 3 through 7 of this review were obtained with the determinantal 
QMC and the DMRG techniques. Here, these two approaches to the many-body problem will be 
described briefly. 

10.1 Determinantal QMC technique 

The determinantal QMC algorithm used in obtaining the QMC data shown here is described in 
Ref. [White et al. 1989b]. Reviews of this technique can be found in Refs. [Scalapino 1993, 

Muramatsu 1999]. The basic idea of this approach is due to Blankenbecler, Scalapino and Sugar 
[Blankenbecler et al. 1981]. Here, the purpose is to calculate the expectation value of an operator 
O at finite temperature T in the grand canonical ensemble, 

(0) = ^TrOe-^^ (10.1) 

where 

Z = Tre-^^. (10.2) 
It is convenient to take the Hubbard hamiltonian zs, R = K -\-V where 

^ = + ^-c-) - Yj^ni^ -I- ni{) (10.3) 

(*i>,o- » 

and 

i 

After discretizing the inverse temperature into L time slices, /? 
is used to rewrite the partition function as 

Z = Tre-^^^" ^ Tr {e'^^^ e'^^^)^ . (10.5) 



(10.4) 

= LAt, the Trotter approximation 
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Next, the interaction term V, which is quartic in the fermion operators, is reduced to a bihnear 
form with the Hubbard-Stratonovich transformation 



„-AtU/4: 

-Arl/K|-l/2)(n^-l/2) ^ £ 



g-ArS'i£A(ni|-n,|) (10.6) 



where cosh(AAr) = exp({7AT/2). Here, at each lattice site i and for each time shce I, an Ising 
spin Sii = ±1 is used for the decouphng. This form of the Hubbard-Stratonovich decouphng was 
introduced by Hirsch [Hirsch 1985]. After integrating out the fermion degrees of freedom, the 
partition function is given by 

Z=Y, detMmSi,})detM^{{Sa}) (10.7) 

{Su} 

where 

M'^iiSie}) = 1 + BlBl_,...B^ (10.8) 

and 

= e-'^^^'^We-^^'^ (10.9) 

with electron spin a = ±1. Here, the summation is over all configurations {Su} of the Ising fields, 
/ is the N X N unit matrix, B^'s are N x N matrices where = SijSa and k is the matrix 

representation of the kinetic energy operator, Eq. (10.3). Similar expressions can be obtained for 



the expectation value (O). For instance, the equal-time single-particle correlation function (cjo-cj,^) 
can be expressed as 

(ci<x4,J = I E ^-(J'-?"! {S^i})detM^{Sa})detM^i{Si,}) (10.10) 

{Su} 

where 

GAj,j'■,{S^e}) = [{I + BlBl_^...Bl)-^].^,. (10.11) 
The summation over the Ising spin variables in Eq. ( 10.10| ) is evaluated using Monte Carlo 



sampling techniques. In this way, Eq. ( |10.10 ) becomes 



(c,<x4j = (G^(j,/;{5,a))p (10.12) 

where (...)p is evaluated by averaging over spin configurations {Si^} generated with the probability 
distribution 

P{{S'U]) = |detMT({Si,})detMH{5,a). (10.13) 

Similarly, the time-dependent single-particle Green's function {cja{T)cjia{T')) is calculated by av- 
eraging Ga{j,T;j' ,T';{Si£}), which can also be expressed in terms of -B/'s. 

The Hubbard-Stratonovich transformation maps the interacting fermion problem on the 2D 
lattice to that of non-interacting fermions in a random Ising field on a 3D lattice, the third dimension 
being the imaginary time axis r. At this stage, for a given spin configuration {Su}, the Wick's 
theorem applies to the higher order correlation functions. Hence, the higher-order correlation 
functions can be obtained by averaging over the products of the single-particle propagators in 
random Ising fields. For instance, the two-particle correlation function 

(rc|(z4,r4)q(«3,r3)c|(i2,r2)c|(ii,ri)), (10.14) 
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which is used in extracting the particle-particle reducible interaction T in Section 6, can be obtained 
from 

(G-f (u, T4; ii, Ti; {Sii})Gi (is, T3; ^2, r2; {Sie}))p. (10.15) 

For (n) = 1, the product of the fermion determinants, det {{Su}) det M-''({5'j£}), is positive. 
However, away from half- filling and for U > 0, this product is not positive for all configurations 
{Sii}. In this case, one can use the probability distribution 

p.r^ .^ ^ \detMH{S,e})detM^{{S,e})\ 

''^^ E{5«} |detMT({S,a)detMH{S,a)|' ^'''''^ 

in order to calculate the expectation value (O) with 

(0({M)sign({g,a))p .,.,7. 
^"^^ (sign({S,4))p ■ ^ ' ' 

Here, sign({5j£}) is the sign of detAfl^dS'j^}) det M-'^dS'j^}) and the average {■■■)p is calculated with 
the probability distribution P. The denominator in Eq. ( |10.17 ), 



(sign) = (sign({5,a))p, (10.18) 

becomes small when the number of configurations with positive signs is close to that with the 
negative signs. It has been shown that away from half-filling (sign) decreases exponentially as T 
decreases [Loh et al. 1990]. This causes large statistical errors in (O) and it is the cause of the 
sign problem. The temperature and the doping regime of the 2D Hubbard model which cannot be 
studied with QMC because of the sign problem is shown at the end of Section 6. In spite of the sign 
problem, useful information has been extracted from the 2D and the 2-leg Hubbard models with 
QMC simulations. Currently, the search for ways of removing the sign problem is an important 
field of study. Improvements in the sign problem which decreases the temperatures accessible to 
QMC simulations by even a factor of two would be valuable. 



10.2 DMRG technique 

The DMRG technique which was invented by White represents a significant development in the 
application of the renormalization group ideas to interacting lattice models [White 1992 and 1993]. 
Currently, it is the numerical method of choice for studying the ground state properties of quasi- 
ID interacting systems. The brief discussion of this technique given here follows closely that of 
Ref. [White 1993]. For simplicity, the following discussion will be for a ID lattice. 

One of the differences between the DMRG and the standard real-space RG approach is in the 
treatment of the boundary conditions in the basic blocking procedure. In the standard approach for 
a finite ID system, the chain is broken into finite blocks Bi of i sites. Here, first the Hamiltonian 
for two identical blocks is diagonalized and then the lowest eigenstates are kept to form a new 
approximate Hamiltonian describing a larger block i?2£ with 2i sites. This procedure is illustrated 
in Fig. 10.1(a). White and Noack have shown that this was of blocking introduces large errors for 
a model of a free particle on a ID lattice because of the way the boundary conditions are treated 
during the blocking [White and Noack 1992]. In the standard RG scheme, neglecting the connection 
of the two blocks to the neighboring blocks corresponds to setting the wave function of the particle 
to outside of the blocks. Consequently, the low lying states from the previous iteration cause a 
'kink' in the wave function in the middle of the enlarged block ^2^. Hence, in order to accurately 
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(a) (b) 

Figure 10.1: (a) Standard blocking scheme used for real-space RG for a ID system. Here, the 
Hamiltonian for two blocks Bg each with £ sites is diagonalized and truncated to form an approx- 
imate Hamiltonian for the new block Bi^^ with 2i sites, (b) DMRG blocking scheme for a finite 
ID system using a superblock which is composed of blocks B£ and S^, and two additional sites in 
between. Here, first the Hamiltonian for the superblock is diagonalized and the density matrix is 
formed for the enlarged block B'^_^^. The Hamiltonian for block B'^_^_^ is then expressed in a reduced 
basis composed of the leading eigenstates of this density matrix. In the next iteration, the block 
B'^_^_-^ replaces B^. 



represent the states of block i?2£, it is necessary to use almost all of the states of block Bi, and any 
truncation of these states introduces large errors in the RG process. White and Noack have noted 
that using a combination of boundary conditions during the blocking fixes this deficiency for a free 
particle on a ID lattice. They have also shown that an alternative approach is to use a superblock 
which is composed of more than two blocks. In this way, two of the blocks in the superblock are 
used to form a larger block for the next iteration while the effect of the other blocks is to apply a 
variety of boundary conditions. The DMRG technique uses a superblock during the RG iteration, 
which is illustrated in Fig. 10.1(b) for a ID lattice of L sites. This superblock is composed of two 
blocks with i and i' = L — i — 2 sites and two additional sites in the middle, which can be considered 
as additional blocks. The block B^ and its neighboring lattice site are called the 'system', while 
the rest of the superblock is called the 'environment'. With each iteration, the system block gets 
enlarged by one lattice site. 

Another novel feature of the DMRG technique is the use of the density matrices to choose the 
states which are to be kept during the iteration [White 1992 and 1993]. In the standard RG ap- 
proach, the lowest m eigenstates of the Hamiltonian for two blocks is used in forming the truncated 
Hamiltonian for the larger block, as illustrated in Fig. 10.1(a). This would be a reasonable approx- 
imation for a model where the coupling between the blocks is weak. However, for an interacting 
system such as the Hubbard model, each block is strongly coupled to its environment. In this 
case, White has shown that it is better to use during the truncation the eigenstates of the density 
matrix of the system, and, in fact, the optimal states to be kept are the eigenvectors of the density 
matrix of the system with the largest eigenvalues. The use of the superblock and the density matrix 
formulation are the important new ideas used by the DMRG technique. 
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At this point, the DMRG algorithm for a finite ID lattice can be summarized as follows: (1) The 
first step is to diagonalize the Hamiltonian for the superblock, which is illustrated in Fig. 10.1(b), 
using exact diagonalization techniques in order to find the ground state ipij . Here, the index i refers 
to the system part of the lattice and j refers to the environment. (2) Next, the reduced density 
matrix for block -B^ , ^ is calculated by tracing over the environment. 



and then p is diagonalized. (3) At this stage, the Hamiltonian for B'^,^-^ can be transformed to a 
reduced basis which consists of the m leading eigenstates of p. Typically a few hundred eigenstates 
of p are kept. (4) In the next iteration, this approximate Hamiltonian for B'^^^ is used in place of 
the Hamiltonian for Bg. 

With this procedure, the system block is built up by adding one lattice site at a time as the 
iteration continues form one end of the lattice to the other end. With each sweep through the 
lattice, a better approximation is obtained for the Hamiltonian of the block B^. In addition, 
during the sweeps, the Hamiltonians of the system blocks from the previous sweep are used as the 
Hamiltonians for the environment blocks. For ID lattices, reflections of i?^ can be used in place 
of B^ in order to build up the superblock to the proper size during the first sweep. In order to 
extend this algorithm to 2D, a row of neighboring sites can be added to the system block at each 
iteration, or a single site can be added using a connected ID path through the 2D lattice. For 
higher dimensional lattices, an empty environment block B^ is used in the first sweep. 

The DMRG technique is more accurate when the number of connections between the system 
and the environment blocks is minimized. Hence, usually open boundary conditions are employed. 
When periodic boundary conditions are used, the accuracy decreases. For 2D lattices, the accuracy 
of DMRG also decreases as a function of the width of the lattice. However, if a larger number of the 
eigenstates of p are kept during the truncation, the accuracy of the algorithm increases. Finally, 
the DMRG is most accurate for large U and least accurate for = 0. 

With this technique, various equal-time correlation functions for the 2-leg Hubbard ladder have 
been calculated on lattices up to 2 x 32 in size [Noack et al. 1996]. For the 3-leg Hubbard ladder 
there are DMRG data for a 3 x 12 lattice [Bonca et al. 2000]. For the t-J model, larger systems can 
be studied with this technique. In this case, 2D clusters up to 28 x 8 in size have been studied using 
cylindrical boundary conditions [White and Scalapino 1998a]. Currently, the DMRG technique is 
used in a wide variety of fields, and the efforts to develop new algorithms to study larger clusters 
continues. In addition, there is continuing work to extend the DMRG technique in order to calculate 
the dynamical correlation functions. 




(10.19) 
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